# Non-Empirical Interactions for the Nuclear Shell Model: An Update

S. Ragnar Stroberg1,2,3 Heiko Hergert4 Scott K. Bogner4 and Jason D. Holt1
1TRIUMF, Vancouver BC, Canada V6T 2A3  2Physics Department, Reed College,
Portland, OR 97202, USA  3Department of Physics, University of Washington,
Seattle WA, USA  4 Facility for Rare Isotope Beams and Department of Physics &
Astronomy, Michigan State University, East Lansing MI 48824, USA

###### Abstract

The nuclear shell model has been perhaps the most important conceptual and
computational paradigm for the understanding of the structure of atomic
nuclei. While the shell model has been predominantly used in a
phenomenological context, there have been efforts stretching back over a half
century to derive shell model parameters based on a realistic interaction
between nucleons. More recently, several ab initio many-body methods—in
particular many-body perturbation theory, the no-core shell model, the in-
medium similarity renormalization group, and coupled cluster theory—have
developed the capability to provide effective shell model Hamiltonians. We
provide an update on the status of these methods and investigate the
connections between them and potential strengths and weaknesses, with a
particular focus on the in-medium similarity renormalization group approach.
Three-body forces are demonstrated to be an important ingredient in
understanding the modifications needed in phenomenological treatments. We then
review some applications of these methods to comparisons with recent
experimental measurements, and conclude with some remaining challenges in ab
initio shell model theory.

###### Contents

  1. 1 Introduction
    1. 1.1 The (Long) Road Towards _Ab Initio_ Shell Model Interactions
    2. 1.2 Fresh Perspectives
      1. 1.2.1 Chiral Effective Field Theory
      2. 1.2.2 The Renormalization Group
    3. 1.3 Current Status of the _ab initio_ Shell Model
    4. 1.4 Organization of this Work
  2. 2 Microscopic effective interactions
    1. 2.1 Quasidegenerate perturbation theory and the $\hat{Q}$-box resummation
      1. 2.1.1 Extended Kuo-Krenciglowa approach
    2. 2.2 Okubo-Lee-Suzuki transformation
    3. 2.3 In-medium similarity renormalization group
      1. 2.3.1 Normal ordering
      2. 2.3.2 Magnus formulation
    4. 2.4 Shell model coupled cluster
  3. 3 Comparison of various approaches to effective interactions
    1. 3.1 Formal effective interaction theory
      1. 3.1.1 General considerations
      2. 3.1.2 Lee-Suzuki
      3. 3.1.3 Generalized Bloch equation
      4. 3.1.4 Shell model coupled cluster
      5. 3.1.5 Canonical perturbation theory
      6. 3.1.6 Fock-space canonical perturbation theory
      7. 3.1.7 Unitary coupled cluster
      8. 3.1.8 Magnus IMSRG
    2. 3.2 Approximation schemes
      1. 3.2.1 Fock-space cluster truncation
      2. 3.2.2 Valence cluster expansion
  4. 4 Three-body forces and the connection with phenomenological adjustments
    1. 4.1 Ensemble normal ordering
    2. 4.2 Mass dependence of the effective interaction
  5. 5 Applications
    1. 5.1 Ground and Excited States of $sd-$Shell Nuclei
    2. 5.2 The Calcium Region
    3. 5.3 Heavy Nickel and Light Tin
  6. 6 Current challenges
    1. 6.1 Electromagnetic transitions
    2. 6.2 The Intruder-state Problem
      1. 6.2.1 The intruder-state problem in perturbation theory
      2. 6.2.2 The intruder-state problem in IMSRG
      3. 6.2.3 Toy model for the intruder-state problem
  7. 7 Other developments
    1. 7.1 EFT for the shell model
    2. 7.2 Uncertainty quantification
    3. 7.3 Coupling to the Continuum
  8. 8 Concluding remarks
  9. A IMSRG flow equations
  10. B Canonical perturbation theory to second order
  11. C Integration of Magnus flow equation to second order

##  1 Introduction

Since its introduction by Goeppert-Mayer and Jensen almost 70 years ago [1, 2,
3], the nuclear shell model has provided the primary conceptual framework for
the understanding of nuclear structure. Its central idea is that protons and
neutrons inside a nucleus move within a self-consistently generated mean
field. This mean field can be approximated by Woods-Saxon or harmonic
oscillator potentials, provided a strong spin-orbit component is added. With
the inclusion of the latter, the spectrum of single-particle energies exhibits
a shell structure that explains the experimentally observed magic proton and
neutron numbers. In this picture, the low-lying structure of most nuclei
results from the interactions between configurations of a few _valence_
particles on top of an inert core via a residual nuclear force.

From the beginning, it was hoped that the shell model and its residual force
could be derived from basic building blocks, in particular the interaction
between free protons and neutrons. Despite early successes, this proved to be
an enormously difficult task [4, 5, 6, 7, 8]. Meanwhile, increasingly
elaborate empirical interactions were used with spectacular success to
describe a vast array of experimental data [9, 10, 11]. In modern language,
this is a clear signal that the shell model contains the relevant degrees of
freedom to describe most (if perhaps not all) phenomena observed in low-energy
nuclear structure.

In the last two decades, a combination of computational and theoretical
advances have provided fresh perspectives and opportunities for the systematic
derivation of shell model interactions from realistic nuclear forces111By
“realistic”, we mean interactions that are rooted to some extent in Quantum
Chromodynamics and accurately describe few-body scattering and bound-state
data., without the need for ad hoc phenomenological adjustments. The primary
aim of this review is to describe these developments. The story of the shell
model, and of microscopically-derived effective interactions in particular, is
long and full of false starts, dead ends, accidental successes, circling back,
and rediscovery of old wisdom in new language. We are not equal to the task of
writing an authoritative history, and we have not attempted to do so. Indeed,
we barely touch on recent developments in empirical shell model interactions
or the computational aspects of configuration interaction calculations.
Readers interested in such techniques should consult the existing texts on the
subject [12, 13, 14, 15]. For the status of phenomenological approaches, we
refer to several thorough reviews [16, 17, 18, 19, 20, 21, 22]. For more
details on effective interaction theory, we recommend several texts [23, 24,
25], and reviews [26, 27, 28, 17, 29, 30].

###  1.1 The (Long) Road Towards _Ab Initio_ Shell Model Interactions

Soon after the publication of the seminal works by Goeppert-Mayer and Jensen,
first parameterizations of the residual nuclear force based on fits to
experimental spectra began to appear (see, e.g., [31, 32, 33, 34, 9]).
Dissatisfied with such approaches because they hide the link between the
effective valence-space interaction and the underlying nucleon-nucleon
interaction, Brown, Kuo and collaborators pioneered the program to
systematically derive and explore this connection [5, 6, 7]. In the 1960s,
considerable progress was made in the construction of realistic nucleon-
nucleon interactions that described NN scattering phase shifts with high
quality [35, 36]. Following a similar strategy as theoretical studies of
nuclear matter, Kuo and Brown used Brueckner’s $G$ matrix formalism [37, 38,
39] to treat the strong short-range correlations induced by these forces, and
included up to second-order diagrams in $G$ to account for important core-
polarization effects [5]. Their work culminated in the publication of
Hamiltonians for the $sd$ and $pf$ shells [7, 8]. While first applications
were successful, Barrett, Kirson, and others soon demonstrated a lack of
convergence in powers of the $G$ matrix [40, 41], and more sophisticated
treatments with RPA phonons and non-perturbative vertex corrections destroyed
the good agreement with experiment [42, 43, 44]. Moreover, Vary, Sauer and
Wong found that the use of larger model spaces for intermediate-state
summations in diagram evaluations also degraded the agreement with
experimental data [45]. Around the same time, Schucan and Weidenmüller [46,
47] demonstrated that the presence of low-lying states of predominantly non-
valence character (“intruder” states) render the perturbative expansion for
the effective interaction divergent.

Because of these developments, enthusiasm for the perturbative approach to
deriving the effective interaction dwindled [48]. While Kuo and collaborators
pursued the $\hat{Q}$-box (or folded-diagram) resummation of the perturbative
series [49, 26, 50, 24, 28], the majority of efforts in shell model theory
were instead focused on the construction and refinement of empirical
interactions [51, 52, 53, 10, 54, 55]. Large-scale calculations with such
interactions yielded impressive agreement with available nuclear data, and
even provided predictive power. An example is the “gold standard” universal
$sd$ shell interaction (USD) [10, 56, 57], which achieves a root-mean-square
deviation from experimental levels of merely 130 keV throughout the $sd$
shell.

Meanwhile, serious efforts were undertaken to develop approaches that
circumvent the problems plaguing the effective interaction methods by starting
from the “bare” nuclear interactions and treating all nucleons as active
particles. Prominent examples are coordinate-space Quantum Monte Carlo (QMC)
techniques [58, 59, 60], or the no-core shell model (NCSM) [61, 62, 63]. The
late 1970s also saw a wave of nuclear Coupled Cluster calculations [64, 65,
66], which use systematic truncations to solve the Schrödinger equation at
polynomial cost, as opposed to the exponential scaling of the NCSM. By the
mid-1990s, computational advances made quasi-exact calculations for nuclei
feasible [67, 68, 61, 62, 58]. However, such calculations were limited to
light nuclei by the sheer numerical cost of coordinate-space QMC, and the slow
convergence of configuration-space methods with realistic nuclear
interactions.

###  1.2 Fresh Perspectives

Since the turn of the millennium, nuclear theory has undergone an important
philosophical shift with the the adoption of renormalization group (RG) and
effective field theory (EFT) concepts. These tools provide a systematic
framework for exploring long-existing questions pertaining to the phase-shift
equivalency of significantly different nucleon-nucleon interactions or the
origin and importance of three-nucleon forces (see, e.g., [69, 70]). Most
noteworthy for the present work is the clarification of the issues that led to
the failure of the aforementioned $G$-matrix based approaches, and the
capability to reconcile the shell model, which is based on an (almost)
independent-particle picture analogous to that of atomic physics, with the
notion that strong correlations are induced by (most) realistic NN
interactions.

####  1.2.1 Chiral Effective Field Theory

The essential idea of chiral effective field theory (EFT) for the nuclear
force is that processes relevant for nuclear structure do not resolve the
details of short-range interactions between nucleons222“Short” in this context
refers to distances $r$ where $k_{F}r\lesssim 1$, with $k_{F}\approx
1.4\,\text{fm}^{-1}$ the Fermi momentum at saturation density.. There are in
fact infinitely-many different potentials, differing at short distances, which
all describe low-energy observables equally well. This is good news, because
we can take advantage of this arbitrariness and parameterize the short-range
physics in a convenient way, e.g., through a series of contact interactions.
At long distances, the approximate chiral symmetry that chiral EFT inherits
from Quantum Chromodynamics (QCD) dictates that interactions are described by
(multi-)pion exchange.

In a pioneering work [71, 72], Weinberg developed effective Lagrangians to
model the interaction between nucleons in terms of pion exchange and contact
interactions, with increasingly complicated contributions suppressed by powers
of a nucleon’s typical momenta $Q\sim k_{F}$ or the pion mass $m_{\pi}$ over
the breakdown scale $\Lambda_{\chi}$ of the EFT, $(Q/\Lambda_{\chi})^{n}$.
This provided a framework to treat two-body forces consistently with three-
and higher-body forces, as well as a natural explanation for the relative
importance of these terms (see e.g. [73, 74, 75] for recent reviews).
Moreover, nuclear transition operators can be derived in a consistent fashion
by coupling the chiral Lagrangian to the electroweak fields (see, e.g., [76,
77, 78, 79, 80]). Despite a number of subtle issues which persist to this day,
several families of chiral two- plus three-nucleon interactions have been
developed [81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92] that reproduce low-
energy observables with an accuracy comparable with phenomenological
potentials. These interactions have become the standard input for modern
nuclear theory.

####  1.2.2 The Renormalization Group

The renormalization group (RG), in particular in the formulation developed by
Wilson [93, 94, 95], is a natural companion to any effective field theory. As
discussed above, an EFT requires a cutoff $\Lambda$ which delineates between
“resolved” and “unresolved” physics. The specific form and location of the
cutoff (the “scheme” and “scale”) is arbitrary and observables for momenta
$Q\ll\Lambda$ should not depend on this choice. Consequently, there are an
infinite number of equivalent theories which differ only in scheme and scale.
The RG smoothly connects such equivalent theories.

RG methods debuted in low-energy nuclear physics around the turn of the
millenium [96, 97, 98, 99, 100, 101, 102, 103], finally providing a systematic
framework which formalized ideas that had been discussed in the nuclear
structure community since the 1950s. For instance, both hard- and soft-core NN
potentials can be devised which reproduce NN scattering data, but nuclear
matter calculations found that soft potentials do not produce empirical
saturation properties and so soft potentials were disfavored [70]. The missing
piece in the saturation puzzle was the connection between the off-shell NN
interaction and the 3N interaction, as demonstrated formally by Polyzou and
Glöckle [104]. Of course, hard-core potentials are much more difficult to
handle in many-body calculations, necessitating the use of Brueckner’s
$G$-matrix [38, 105, 106, 107] to deal with correlations due to the short-
range repulsion.

From the RG perspective, the hard- and soft-core potentials are related by an
RG transformation which leaves NN scattering observables unchanged, but shifts
strength into _induced_ 3N (and higher) interactions. Neglecting these induced
terms means that observables involving more than two particles will no longer
be preserved. This mechanism provides an explanation of the Phillips line
[108] and the Tjon line [109], which describe correlations between few-body
observables calculated using different phase-equivalent NN interactions.

In the context of this work, the RG provides a simple explanation of the
observation that — after being processed by the Brueckner $G$-matrix machinery
— various NN potentials produce very similar spectroscopy [110, 111], as long
as they reproduce NN scattering data. This can be understood as an indication
that the $G$-matrix approximately “integrates out” the short-distance physics
of the different potentials, leaving the universal long-distance physics.
However, Bogner _et al._ have shown that the $G$-matrix can retain significant
coupling between off-shell low- and high-momentum modes, rendering it non-
perturbative [102]. This explains why the historical efforts to construct the
effective interaction perturbatively from $G$-matrices were bound to fail. In
contrast, methods like the Similarity RG (SRG) [112, 113, 102] achieve a more
complete decoupling of the short-distance physics and render the resulting
transformed NN+3N interaction suitable for perturbative expansions [114, 102,
115, 116]. The SRG has become the tool of choice in nuclear theory for
decoupling low and high momenta because it also provides straightforward means
to track induced many-body forces [117, 118, 119, 120], and construct
consistently transformed observables [121, 122, 123].

###  1.3 Current Status of the _ab initio_ Shell Model

Over the past decade, the adaptation of EFT and SRG methods has greatly
extended the reach of _ab initio_ nuclear many-body theory across the nuclear
chart. Simply put, the most convenient scale for formulating a theory of
nuclear interactions is often not the most convenient scale for solving that
theory. The SRG connects one scale to the other and greatly improves the
convergence behavior of nuclear many-body calculations in the process. Large-
scale diagonalization methods like the NCSM can be used in the lower
$sd-$shell [63, 124, 125, 120], and systematically truncated methods like
Self-Consistent Green’s Functions, Coupled Cluster (CC) and the In-Medium SRG
(IMSRG) can even be applied to nuclei as heavy as tin [126, 127, 128, 129].
While SRG-evolved interactions cannot be used easily in QMC due to their
nonlocality, new families of local chiral interactions yield encouraging
results in such applications [87, 88, 90, 60].

Soon after their introduction to nuclear physics, EFT and RG methods also
revitalized efforts to systematically derive shell model interactions [130,
131]. From a practical perspective, this offered a convenient way to confront
RG-evolved chiral two- plus three-nucleon interactions with the wealth of
available spectroscopic data, using existing shell model codes [132, 133, 134,
135, 136, 137, 138, 139, 140, 141, 142, 143, 144]. At the conceptual level,
these interactions validate the independent-particle picture underlying the
shell model. They provide sufficient binding already at the mean-field level,
and allow us to use it as the starting point for the treatment of
correlations, either through rapidly converging non-perturbative many-body
methods [145, 146, 127], or possibly even through finite-order perturbation
theory [116, 147]. Furthermore, novel approaches like the valence-space IMSRG
(VS-IMSRG) or shell model Coupled Cluster (SMCC), both discussed below,
provide both the conceptual framework and practical tools to relate no-core
and valence-space methods, as shown by the consistent ground- and excited-
state results obtained thus far (see [143, 129] and Sec. 5). Thus, the end of
the long and winding road to _ab initio_ shell model interactions appears to
be in sight, though challenges remain, see Secs. 6 and 8.

###  1.4 Organization of this Work

This work is organized as follows: In Sec. 2, we introduce common approaches
to the construction of shell model interactions, from the traditional many-
body perturbation theory and the Okubo-Lee-Suzuki method to the valence-space
IMSRG and shell model CC. In Sec. 3, these approaches are compared within a
common formalism to illuminate the relations between them. Section 4 discusses
the role of three-nucleon forces in the shell model context, and relates
modifications of (semi-)empirical interactions that are supposed to capture
such effects to the more systematic treatment of these forces in modern
approaches. In Sec. 5, we highlight selected applications of _ab initio_ shell
model interactions. Section 6 describes the main challenges we are facing
today, and analyzes them primarily from the perspective of the VS-IMSRG. New
developments like a direct EFT expansion for shell model interactions and a
novel uncertainty quantification effort are touched upon in Sec. 7. Section 8
provides concluding remarks alongside a list of take-away messages that
summarize the key aspects of modern _ab initio_ shell model calculations, and
clarify common misconceptions. Certain technical details are collected in the
appendices.

##  2 Microscopic effective interactions

The general problem of effective interaction theory is the following: Given a
Hamiltonian $H$ expressed in a large (typically intractable) Hilbert space
$\mathscr{H}$, we wish to obtain an effective Hamiltonian $H_{\mathrm{eff}}$
which acts in a smaller (tractable) Hilbert space
$\mathscr{H}_{\mathrm{model}}$, but reproduces a subset of the eigenstates of
the large Hilbert space.

$\underset{\textrm{Full space Schr\"{o}dinger
eq.}}{\underbrace{H|\Psi_{n}\rangle=E_{n}|\Psi_{n}\rangle}}\hskip
10.00002pt\Rightarrow\hskip 10.00002pt\underset{\textrm{Model space
Schr\"{o}dinger
eq.}}{\underbrace{H_{\text{eff}}|\psi_{n}\rangle=E_{n}|\psi_{n}\rangle}}$ (1)

In the context of the nuclear shell model, the large Hilbert space will
consist of Slater determinants of single-particle states, typically harmonic
oscillator eigenstates. The number of single-particle states included should
be sufficient to obtain convergence. The smaller Hilbert space
$\mathscr{H}_{\mathrm{model}}$ is defined by splitting the single-particle
states into three categories—core, valence, and excluded333In chemistry,
valence orbits are usually called active states, while excluded orbits are
referred to as virtual states.— taking the subset of Slater determinants for
which all core orbits are occupied, and all excluded orbits are unoccupied.

We note in passing that different partitionings of the Hilbert space can be
used for other purposes. For example, choosing $\mathscr{H}_{\mathrm{model}}$
to be a one-dimensional space corresponds to the single-reference many-body
methods for treating the ground state of closed-shell nuclei. Alternatively,
defining $\mathscr{H}_{\mathrm{model}}$ in terms of low-momentum states leads
to methods for “softening” an interaction, such as $V_{\mathrm{low}k}$ [99] or
the similarity renormalization group (SRG) [100].

An important practical requirement on $H_{\mathrm{eff}}$ is that it should
obey a rapidly converging cluster expansion, schematically

$|V_{2N}|\gg|V_{3N}|\gg|V_{4N}|\ldots$ (2)

where $V_{2N}$ is the two-body potential, $V_{3N}$ is the three-body
potential, etc. and the vertical bars indicate some measure of size or
importance. This property is essential to the feasibility of large-scale shell
model diagonalizations. Nowadays, such calculations can handle basis
dimensions upwards of $10^{9}$ Slater determinants [18, 20, 148], which would
be impossible if the full matrix needed to be stored. The limitation to two-
body (or possibly three-body) interactions yields a sparse matrix which can be
treated efficiently by, e.g., the Lanczos or Davidson methods [149, 150, 151].

It is worth taking a moment here to motivate why we would expect the cluster
expansion to be valid in nuclei, and to consider where it might run into
trouble. As mentioned in section 1.2.1, chiral effective field theory
naturally generates a hierarchy of the type shown in Eq. (2), with many-body
interactions suppressed by increasing powers of the ratio of low to high
scales. On the other hand, for a system of $A$ nucleons, the importance of an
$n$-body term grows combinatorially, accounting for all the different
combinations of $n$ particles which can interact, suggesting that for heavy
nuclei, many-body forces will dominate. For large $A$, this grows as
$A!/n!(A-n)!\sim A^{n}$.

Fortunately, we are saved by the short range of the nuclear force and the
relatively low saturation density of nuclear matter [152]. Roughly, each
nucleon does not interact with all other nucleons, but instead only with the
other nucleons within some interaction volume
${\mathcal{V}\sim\tfrac{4\pi}{3}r_{\mathrm{int}}^{3}}$ with $r_{\mathrm{int}}$
the range of the interaction. At density $\rho$, the expectation value of an
$n$-body force will scale as $\langle
V_{n}\rangle\sim(\rho\mathcal{V})^{n-1}$. With $\rho\lesssim
0.16\,\text{fm}^{-3}$, the cluster hierarchy will be maintained so long as
$r_{\mathrm{int}}\lesssim 1\,\text{fm}$. This is satisfied for the short-range
terms in a chiral force, which are characterized by a cutoff scale on the
order of 0.5 fm or less. The long pion-exchange tail, with range $r_{\pi}\sim
1.5\,\text{fm}$, is not obviously suppressed or enhanced by the density,
though its contribution to bulk binding is somewhat suppressed due to its
spin-isopin structure averaging to zero in symmetric spin-saturated matter.
However, if we use Wick’s theorem to express operators in normal ordered form
with respect to a finite-density reference (see sections 2.3.1 and 4.1), then
the appropriate density is not saturation density $\rho$, but the
_quasiparticle_ density which will typically be significantly smaller [152].

When we derive an effective interaction, we are eliminating degrees of
freedom, namely orbits outside of the valence space. For high-lying orbits,
the relevant interaction matrix element will be dominated by high-momentum
(short-distance) components, in which case the above argument holds and the
induced terms should still exhibit a cluster hierarchy [153]. However, for
excitations near the fermi surface we have no short-distance argument. Indeed,
as we discuss in section 6.1, low-lying collective excitations can be a source
of trouble.

It is important to keep in mind that the above argument only holds if the
observable in question can be expressed in terms of connected diagrams (i.e.
one can trace a continuous path through the diagram between any two points on
it). If a four-body term consists of two disconnected two-body terms, then
there is no reason why all four particles would need to be within some
interaction volume. We return to this point in Sec. 3.2.2.

In the remainder of this section, we describe the most popular approaches to
deriving effective interactions for the nuclear shell model, using the
notation that appears in the nuclear physics literature. In section 3, we
treat these methods within a more general framework to illuminate the
relationships between them.

###  2.1 Quasidegenerate perturbation theory and the $\hat{Q}$-box resummation

Let us introduce the projection operator $P$ and its complement $Q$ such that
$P\mathscr{H}P=\mathscr{H}_{\mathrm{model}}$ and $P+Q=1$. If we assume
$|\psi_{n}\rangle=P|\Psi_{n}\rangle$, that is, that the eigenstates of the
effective Hamiltonian are just the projection of the full eigenstates onto the
model space, then the effective Hamiltonian should satisfy

$PH_{\mathrm{eff}}P|\Psi_{n}\rangle=E_{n}P|\Psi_{n}\rangle\hskip
10.00002pt,\hskip 10.00002ptQH_{\mathrm{eff}}P=0.$ (3)

Straightforward manipulation then yields the Bloch-Horowitz energy-dependent
effective Hamiltonian [154, 155]

$H_{\mathrm{BH}}(E_{n})=PHP+PHQ\frac{1}{E_{n}-QHQ}QHP.$ (4)

One important aspect of Eq. (4) is that the effective interaction depends on
the eigenvalue $E_{n}$, and so it must be solved self-consistently. A second
point is that different valence-space eigenstates will in general not be
orthogonal, because they are eigenstates of different operators. The energy
dependence can be removed by expanding the denominator about some starting
energy $E_{0}$, yielding [156, 138]

$H_{\mathrm{eff}}=H_{\mathrm{BH}}(E_{0})+\sum_{k=1}^{\infty}\frac{1}{k!}\left[\frac{d^{k}}{dE_{0}^{k}}H_{\mathrm{BH}}(E_{0})\right](H_{\mathrm{eff}}-E_{0})^{k}.$
(5)

The expression (5) may also be obtained in the context of time-dependent
perturbation theory [24], or by a similarity transformation combined with an
iterative solution [157] for the decoupling condition (3).

A simplification may be obtained by partitioning the Hamiltonian into a zero-
order piece and a perturbation $H=H_{0}+V$, and assuming that the eigenvalues
of $H_{0}$ in the valence space are degenerate, with energy $\epsilon$. Then
one uses Eq. (5) with $E_{0}=\epsilon$, and
$(H_{\mathrm{eff}}-\epsilon)=V_{\mathrm{eff}}$. In this context, a popular
approach is to define the “$\hat{Q}$-box” [49, 50, 24], indicated
$\hat{Q}(\epsilon)$.

$\hat{Q}(\epsilon)=PVP+PVQ\frac{1}{\epsilon-QHQ}QVP.$ (6)

The operator $\hat{Q}(\epsilon)$ is not to be confused with the projection
operator $Q$. The effective valence space interaction is, in analogy to Eq.
(5),

$V_{\mathrm{eff}}=\hat{Q}(\epsilon)+\sum_{k=1}^{\infty}\frac{1}{k!}\frac{d^{k}\hat{Q}(\epsilon)}{d\epsilon^{k}}(V_{\mathrm{eff}})^{k}.$
(7)

The commonly adopted strategy for evaluating Eq. (7) is to expand the inverse
operator in Eq. (6) perturbatively (typically to second or third order in $V$)
and to solve Eq. (7) self-consistently by iteration, evaluating the
derivatives numerically by calculating $\hat{Q}(\epsilon)$ for several
starting energies $\epsilon$.

A technical point arises because the $\hat{Q}$-box contains one-body pieces,
arising from e.g. core polarization diagrams [5, 6, 7]. For computational
convenience, the one-body part is embedded in the two-body part, with an
accompanying spectator nucleon. This leads to disconnected two-body terms in
Eq. (7) which contain arbitrary numbers of interactions involving one-
particle, but no interactions between the two. These disconnected
contributions can be understood as the dressed one-body part of the effective
interaction embedded into a two-body interaction. They can be removed by
solving Eq. (7) using only the one-body piece of the $\hat{Q}$-box, which is
called the $\hat{S}$-box [158, 159]. The resulting effective one-body
interaction is subtracted off from $V_{\mathrm{eff}}$, leaving only connected
terms. In principle these same diagrams should then be added self-consistently
to the degenerate single-particle energies $\varepsilon_{i}=\langle
i|H_{0}|i\rangle+\hat{S}_{i}(\varepsilon_{i})$. In practice, it is often
discarded and the single particle energies are taken from
experiment444Determining the experimental single-particle energies is
ambiguous, as one must choose which experimental state is “the” shell model
one, or consider an average of multiple states (see, e.g., [10, 56]). As an
exception to this, when three-body forces are included in the normal ordering
approximation (see Secs. 2.3.1 and 4), and the single-particle energies are
computed so that the starting energy $\epsilon$ corresponds to the centroid of
the valence-space single-particle energies, no additional adjustments are
required [136, 137].

####  2.1.1 Extended Kuo-Krenciglowa approach

In general, one would like to be able to use non-degenerate valence orbits,
for example when using a Hartree-Fock basis, or when generating an interaction
for two (or more) major harmonic oscillator shells. Then one should use Eqs.
(4) and (5) rather than Eqs. (6) and (7). This approach has only recently been
explored [160, 161, 138, 144]. In the case of a valence space spanned by two
major oscillator shells, taking harmonic oscillator single-particle energies
can lead to zero energy denominators in Eq. (7), while the parameter $E$ in
Eq. (5) may be chosen to avoid zero denominators. Of course, for a $P$ space
with a very large spread in single-particle energies, one would expect that
the size of certain contributions to $(H_{\mathrm{eff}}-E_{0})^{k}$ in Eq. (5)
would be comparable to the corresponding energy denominators, and the Taylor
series (5) might converge slowly, if at all. However, for a modest spread of
energies ($\sim$10 MeV), this does not appear to be a problem [138, 144].

###  2.2 Okubo-Lee-Suzuki transformation

Another approach to the effective interaction is called the Okubo-Lee-
Suzuki555Perhaps this would more appropriately be called the Okubo-Lee-Suzuki-
Okamoto approach. approach, often employed in conjunction with the NCSM [162,
62, 63], the paradoxical-sounding “NCSM with a core” or Double OLS approach
[163, 164, 165], or more recently, the coupled cluster method [140, 166]
(“coupled cluster effective interaction” (CCEI)). The idea is to obtain a
unitary transformation $\mathcal{U}$ which diagonalizes $H$ in the large
Hilbert space $\mathscr{H}$, so that

$\mathcal{U}H\mathcal{U}^{\dagger}=\mathcal{E}$ (8)

where $\mathcal{E}$ is a diagonal matrix. The effective Hamiltonian
$H_{\mathrm{eff}}$ acting in the smaller Hilbert space
$\mathscr{H}_{\mathrm{model}}=P\mathscr{H}P$ is given by [163]

$H_{\mathrm{eff}}=\frac{\mathcal{U}_{P}^{\dagger}}{\sqrt{\mathcal{U}^{\dagger}_{P}\mathcal{U}_{P}}}\mathcal{E}\frac{\mathcal{U}_{P}}{\sqrt{\mathcal{U}^{\dagger}_{P}\mathcal{U}_{P}}}$
(9)

where $\mathcal{U}_{P}\equiv P\mathcal{U}P$ is the projection of the
transformation $\mathcal{U}$ onto the model space. One can easily confirm that
the transformation in Eq. (9) is unitary, and so the eigenvalues of
$H_{\mathrm{eff}}$ in the model space will be a subset of the eigenvalues of
$H$ in the full space.

So far, this does not appear to be a helpful procedure, because the first step
was to solve the eigenvalue problem in the large Hilbert space, and the goal
of effective interaction theory is to allow applications for which the full
solution is not tractable. The benefit comes when one assumes that the
effective interaction $H_{\mathrm{eff}}$ will also provide a good
approximation for other systems described in the same model space. Assuming a
cluster expansion, one obtains $H_{\mathrm{eff}}$ for a few active particles
in the model space and then applies $H_{\mathrm{eff}}$ to systems with more
active particles.

As a concrete example, consider the system of two particles in the no-core
shell model. This can readily be diagonalized in a large space of many
harmonic oscillator states, say $N\leq N_{\mathrm{max}}=500$, where
$N=2n+\ell$ is the number of oscillator quanta. On the other hand, even a
light nucleus like 6Li with only six particles cannot possibly be diagonalized
in such a space. Using the OLS transformation (9) one can obtain an effective
interaction for a manageable model space (say $N_{\mathrm{max}}=10$) which
exactly reproduces the low-lying eigenvalues of the two-body calculation in
the large space. The application of this effective interaction to 6Li then
provides a reasonable approximation to the eigenvalues one would obtain in the
large space.

Figure 1: A schematic of how the OLS approach obtains the effective
interaction. (a) The original Hamiltonian. (b) The Hamiltonian is diagonalized
by transformation $\mathcal{U}$. (c) The transformation
$\mathcal{U}_{P}^{\dagger}$ approximately inverts $\mathcal{U}$ in the $P$
space and yields the effective interaction $H_{\mathrm{eff}}$.

We note that there are essentially two main assumptions here: _(i)_ the
effective interaction one would obtain if one could apply the OLS procedure
directly to the six-body system has a rapidly convergent cluster expansion
$|V_{2N}|>|V_{3N}|>|V_{4N}|\ldots$, and _(ii)_ the two-body-cluster component
$V_{2N}$ of the full effective interaction $H_{\mathrm{eff}}$ for the six-body
system is well approximated by the effective interaction obtained for the two-
body system. While both of these assumptions are plausible and encouraging
results have been obtained using them, we know of no rigorous proof. Indeed,
as we discuss briefly in Sec. 3.2.2, there is potential cause for concern
related to disconnected diagrams.

###  2.3 In-medium similarity renormalization group

In the in-medium similarity renormalization group (IMSRG)[167, 168, 139, 143,
146], the effective Hamiltonian is also expressed in terms of a unitary
transformation $U$ acting on the initial Hamiltonian

$H_{\mathrm{eff}}=UHU^{\dagger}.$ (10)

In contrast to the OLS approach, the IMSRG transformation is obtained without
solving the eigenvalue problem for a particular many-body system. Instead, it
is parameterized by a continuous flow parameter $s$, and applied to the
Hamiltonian through the flow equation

$\frac{dH(s)}{ds}=[\eta(s),H(s)]\,,$ (11)

where the generator $\eta(s)$ is formally defined as

$\eta(s)\equiv\frac{dU(s)}{ds}U^{\dagger}(s)=-\eta^{\dagger}(s).$ (12)

We split the flowing Hamiltonian $H(s)$ into two pieces, a “diagonal” and an
“off-diagonal” piece

$H(s)=H^{d}(s)+H^{od}(s)$ (13)

such that

$\displaystyle H^{od}(s)$ $\displaystyle=PH(s)Q+QH(s)P.$ (14)

where the projection operators $P$ and $Q$ have the same meaning as in the
previous sections. Our goal is then to devise a generator $\eta(s)$ such that

$\underset{s\rightarrow\infty}{\lim}H^{od}(s)=0$ (15)

and therefore

$\underset{s\rightarrow\infty}{\lim}H^{d}(s)=H_{\mathrm{eff}}.$ (16)

In the language of the renormalization group, $H_{\mathrm{eff}}$ is a fixed
point of the RG flow.

One choice for $\eta(s)$, which is used in the calculations we will describe
here is the White generator [169, 145]

$\eta^{\mathrm{Wh}}(s)\equiv\frac{H^{od}(s)}{\Delta(s)}\,.$ (17)

For present and future use, we have introduced a convenient superoperator
notation (cf. [170]), in which we indicate division of the operator $O$ by a
suitably defined energy denominator $\Delta$ is defined as

$\langle\phi_{i}|\frac{O}{\Delta}|\phi_{j}\rangle\equiv\frac{\langle\phi_{i}|O|\phi_{j}\rangle}{\epsilon_{i}-\epsilon_{j}}$
(18)

which can be thought of as element-wise division. Here
$\epsilon_{i},\epsilon_{j}$ are energies associated with the basis states
$\phi_{i},\phi_{j}$. The quantity $\frac{O}{\Delta}$ itself is an operator
whose Fock-space expression is

$\frac{O}{\Delta}=\sum_{ij}\frac{O_{ij}}{\epsilon_{i}-\epsilon_{j}}a^{\dagger}_{i}a_{j}+\frac{1}{4}\sum_{ijkl}\frac{O_{ijkl}}{\epsilon_{i}+\epsilon_{j}-\epsilon_{k}-\epsilon_{l}}a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}+\ldots$
(19)

Returning to the flow equation, it is clear that if $H^{od}\rightarrow 0$,
then $\eta\rightarrow 0$ and by Eq. (11) we see that
$\frac{dH(s)}{ds}\rightarrow 0$, so $H_{\mathrm{eff}}$ is indeed a fixed point
of the flow. One potential issue with the generator (17) is that a vanishing
energy denominator will cause $\eta$ to diverge. An alternative, also
suggested by White [169] (see also [171]), is

$\eta^{\mathrm{atan}}(s)\equiv\frac{1}{2}\mathrm{atan}\left(\frac{2H^{od}(s)}{\Delta(s)}\right).$
(20)

The arctangent—motivated by the solution of a 2$\times$2 system via Jacobi
rotations—regulates the divergent behavior of Eq. (17) in the presence of
small denominators. The arctangent and division by the energy denominator in
Eq. (20) should be interpreted as operating element-wise, as described above.

Figure 2: A schematic representing of how the IMSRG approach obtains the
effective interaction $H_{\mathrm{eff}}$ by progressively suppressing the off-
diagonal terms of $H$. (a)$s=0$, (b)$s=5$, (c)$s=30$

The IMSRG is formulated in terms of Fock-space operators, and so its
computational cost scales polynomially with the basis size $N$, but not
explicitly with the number of particles being treated. In practical
applications, we truncate all operators at a consistent particle rank to close
the system of flow equations arising from Eq. (11) (see Appendix A). We also
set up the decoupling conditions to be minimally invasive to avoid an
uncontrolled accumulation of truncation errors, as discussed in detail in Ref.
[145]. In VS-IMSRG we perform the decoupling in two stages for this reason,
decoupling the reference state from excitations as in a direct ground-state
calculation before decoupling the valence space in a second evolution (cf.
Secs. 5 and 6).

####  2.3.1 Normal ordering

An important feature of the IMSRG method is the use of operators in _normal-
ordered form_ (see, e.g., [145, 146]). Starting with the free-space
Hamiltonian written as a Fock-space operator with two- and three-body
interactions 666In actual calcualtions, one subtracts off the center-of-mass
kinetic energy, and so the kinetic term has a two-body piece [146]. We neglect
that here for simplicity.

$H=\sum_{ij}t_{ij}a^{\dagger}_{i}a_{j}+\tfrac{1}{4}\sum_{ijkl}V_{ijkl}a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}+\tfrac{1}{36}\sum_{\begin{subarray}{c}ijk\\\
lmn\end{subarray}}V_{ijklmn}a^{\dagger}_{i}a^{\dagger}_{j}a^{\dagger}_{k}a_{n}a_{m}a_{l},$
(21)

we may use Wick’s theorem to express the strings of creation and annihilation
operators in normal order with respect to some reference state $|\Phi\rangle$
[172]. We denote the normal ordering with braces, and the normal order of a
pair of operators is defined so that their expectation in the reference is
zero, i.e.

$\langle\Phi|\\{\\!a^{\dagger}_{i}a_{j}\\}|\Phi\rangle=0.$ (22)

Whether the normal order is $a^{\dagger}_{i}a_{j}$ or $a_{j}a^{\dagger}_{i}$
depends on whether or not the states created and annihilated are present in
the reference $|\Phi\rangle$. If we choose $|\Phi\rangle$ to be a single
Slater determinant such as the Hartree-Fock ground state of the system of
interest, then application of Wick’s theorem allows us to write $H$ as

$H=E_{0}+\sum_{ij}f_{ij}\\{\\!a^{\dagger}_{i}a_{j}\\}+\tfrac{1}{4}\sum_{ijkl}\Gamma_{ijkl}\\{\\!a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}\\}+\tfrac{1}{36}\sum_{\begin{subarray}{c}ijk\\\
lmn\end{subarray}}W_{ijklmn}\\{\\!a^{\dagger}_{i}a^{\dagger}_{j}a^{\dagger}_{k}a_{n}a_{m}a_{l}\\}$
(23)

where the new coefficients can be obtained from the old coefficents by

$\displaystyle
E_{0}=\sum_{a}n_{a}t_{aa}+\tfrac{1}{2}\sum_{ab}n_{a}n_{b}V_{abab}+\tfrac{1}{6}\sum_{abc}n_{a}n_{b}n_{c}V_{abcabc}$
(24) $\displaystyle
f_{ij}=t_{ij}+\sum_{a}n_{a}V_{iaja}+\tfrac{1}{2}\sum_{ab}n_{a}n_{b}V_{iabjab}$
$\displaystyle\Gamma_{ijkl}=V_{ijkl}+\sum_{a}n_{a}V_{ijakla}$ $\displaystyle
W_{ijklmn}=V_{ijklmn}.$

Operators other than the Hamiltonian can be rewritten in the same way.

In Eq. (24), $n_{a}$ is the occupation of orbit $a$ in the reference, i.e.
$n_{a}=\langle\Phi|a^{\dagger}_{a}a_{a}|\Phi\rangle$, and for a Slater
determinant reference, $n_{a}$ is either zero or one. In section 4.1, we
discuss a different choice of reference for which $a$ can have fractional
occupation. One may also use a correlated reference, constructed out of a
linear combination of Slater determinants, in which case one must use the
generalized normal ordering of Kutzelnigg and Mukherjee [173]. This is the
basis of the multireference IMSRG (MR-IMSRG) method that will be used in
ground-state energy comparisons in later sections [125, 146].

The advantage of expressing operators in normal ordered form is that it puts
as much information as possible from the high particle-rank (i.e. many-body)
operators into the lower rank operators. This is evident in Eq. (24), where
the normal-ordered zero-body term $E_{0}$ contains contributions from the free
one-, two- and three-body terms. If the reference $|\Phi\rangle$ is a good
approximation of the exact wave function $|\Psi\rangle$, then the expectation
value
${\langle\Psi|\\{\\!a^{\dagger}a^{\dagger}a^{\dagger}aaa\\}|\Psi\rangle\approx
0}$, and even formally non-vanishing 3N interactions $W_{ijklmn}$ can be
neglected to a good approximation. Consequently, normal ordering may be
thought of as a way to improve the convergence of the cluster expansion
described in the beginning of this section.

####  2.3.2 Magnus formulation

A particularly convenient formulation of the IMSRG approach relies on the
Magnus expansion [174, 175]. The idea is to express the more general unitary
IMSRG transformation as the true exponential of the anti-Hermitian Magnus
operator $\Omega(s)=-\Omega^{\dagger}(s)$. The evolved Hamiltonian can then be
expressed in terms of an infinite series of nested commutators

$\displaystyle H(s)$ $\displaystyle=e^{\Omega(s)}H(0)e^{-\Omega(s)}$ (25)
$\displaystyle=H(0)+[\Omega(s),H(0)]+\frac{1}{2}\left[\Omega(s),[\Omega(s),H(0)]\right]+\ldots$

This formulation of the IMSRG allows for a more transparent comparison with
canonical transformation theory [169, 176] and unitary coupled cluster method
[177] used in quantum chemistry, as well as canonical perturbation theory
[178], where the expansion in Eq. (25) is evaluated perturbatively.

Considering the flow equation (11), we see that under an infinitesimal step
$ds$ we may write

$\displaystyle H(s+ds)$ $\displaystyle=H(s)+[\eta(s),H(s)]ds$ (26)
$\displaystyle=e^{\eta(s)ds}H(s)e^{-\eta(s)ds}$
$\displaystyle=e^{\eta(s)ds}e^{\Omega(s)}H(0)e^{-\Omega(s)}e^{-\eta(s)ds}\,.$

Expressing $H(s+ds)$ in Magnus form as well, we obtain an expression for
$\Omega(s+ds)$:

$e^{\Omega(s+ds)}=e^{\eta(s)ds}e^{\Omega(s)}.$ (27)

Because $\Omega(s)$ and $\eta(s)$ do not in general commute, we use the Baker-
Campbell-Hausdorff formula to take the logarithm on both sides and obtain

$\Omega(s+ds)=\Omega(s)+\eta(s)ds+\tfrac{1}{2}[\eta(s),\Omega(s)]ds+\tfrac{1}{12}[\Omega(s),[\Omega(s),\eta(s)]]ds+\ldots$
(28)

This may be expressed in compact form as (suppressing explicit $s$ dependence)

$\frac{d\Omega}{ds}=\sum_{k=0}^{\infty}\frac{B_{k}}{k!}\mathrm{ad}^{(k)}_{\Omega}(\eta)$
(29)

where $B_{k}$ are the Bernoulli numbers, and the adjoint
$\mathrm{ad}^{(k)}_{\Omega}(\eta)$ signifies a recursively defined nested
commutator:

$\mathrm{ad}^{(k)}_{\Omega}(\eta)=[\Omega,\mathrm{ad}^{(k-1)}_{\Omega}(\eta)]\hskip
10.00002pt,\hskip 20.00003pt\mathrm{ad}^{(0)}_{\Omega}(\eta)=\eta.$ (30)

Fortunately, in most practical applications, only the first few terms of the
infinite series in Eqs. (25) and (28) are important, and so the commutator may
be evaluated iteratively until the size of a given term is below some
numerical threshold (for an example of an exception to this, see Sec. 6.2). A
major practical advantage of the Magnus method is that by solving for
$\Omega(s)$, we can compute arbitrary effective operators besides the
Hamiltonian in a consistent and efficient way (see Sec. 6.1 and Refs. [174,
123]).

###  2.4 Shell model coupled cluster

Figure 3: A schematic depicting how the SMCC approach obtains the effective
interaction $\bar{\bar{\mathcal{H}}}(t)$. Note that, in contrast to the
$H_{\mathrm{eff}}$ of Fig. 2, $\bar{\bar{\mathcal{H}}}$ is non-Hermitian.
(a)$t=0$, (b)$t=5$, (c)$t=30$

First attempts to derive an effective shell model interaction with coupled
cluster methods were similar in spirit to the NCSM-based Double OLS approach
(see Sec. 2.2 and Refs. [162, 163, 164, 165]). Equation-of-motion CC (EOM-CC)
states defined in a space of up to 4p2h excitations are subsequently projected
into the shell model space via the OLS method, yielding the so-called coupled
cluster effective interaction (CCEI) [140, 166]. The cost of the EOM-CC
calculations, however, presented a significant obstacle to widespread
applications of this method. A much more efficient alternative is the recently
introduced shell model coupled cluster (SMCC) method [179], which is
formulated in Fock space and can be viewed as a non-unitary cousin of the
IMSRG.

In the standard coupled cluster method [25, 127], a similarity transformation
is performed to decouple a single closed-shell reference state $|\Phi\rangle$
from all particle-hole excitations

$e^{-T}He^{T}|\Phi\rangle=E_{\mathrm{corr}}e^{T}|\Phi\rangle$ (31)

where $E_{\mathrm{corr}}$ is the correlation energy and $T$ is the cluster
operator which is written as

$T=\sum_{ia}t_{i}^{a}\\{a^{\dagger}_{a}a_{i}\\}+\tfrac{1}{4}\sum_{abij}t_{ij}^{ab}\\{a^{\dagger}_{a}a^{\dagger}_{b}a_{j}a_{i}\\}+\ldots$
(32)

Here the indices $a,b,c\ldots$ denote unoccupied orbits and $i,j,k\ldots$
denote occupied orbits. The similarity-transformed Hamiltonian is written as

$\bar{\mathcal{H}}=e^{-T}He^{T}.$ (33)

Shell model coupled cluster extends the idea by performing a similarity
transformation which decouples a valence space rather than a single
configuration. Denoting this transformation with an $S$, we have

$\bar{\bar{\mathcal{H}}}=e^{-S}\bar{\mathcal{H}}e^{S}\hskip 10.00002pt,\hskip
20.00003ptQ\bar{\bar{\mathcal{H}}}P=0.$ (34)

The operator $S$ is obtained by a flow equation closely mirroring the one used
in IMSRG

$\frac{dS}{dt}=-\eta(\bar{\bar{\mathcal{H}}}(t))$ (35)

where $\eta$ is the generator of the flow. As in the IMSRG, there is
considerable freedom for choosing $\eta$ as long as the decoupling condition
(34) is realized in the limit $t\to\infty$. In Ref. [179], adapted variants of
the White (Eq. (17)) and arctangent generators (Eq. (20)) are used. The
essential difference from the IMSRG formulation is that $T$ and $S$ are not
anti-Hermitian operators, and so the transformation is not unitary, and the
resulting effective Hamiltonian is not Hermitian. The inconvenience of a non-
Hermitian Hamiltonian is compensated by the greater simplicity of the
equations which need to be solved.

##  3 Comparison of various approaches to effective interactions

To investigate how the methods described in section 2 (and a few others) are
related to one another, we consider the general structure of effective
interactions, and show how the above methods sum the perturbation series. For
more details, we refer readers to Refs. [170, 178, 180, 181, 182, 183].

###  3.1 Formal effective interaction theory

####  3.1.1 General considerations

We begin by expressing the effective Hamiltonian in terms of a similarity
transformation of the original Hamiltonian, parametrized as the exponential of
a generator $\mathcal{G}$

$\displaystyle H_{\mathrm{eff}}$
$\displaystyle=e^{\mathcal{G}}He^{-\mathcal{G}}$ (36)
$\displaystyle=H+[\mathcal{G},H]+\tfrac{1}{2}\left[\mathcal{G},[\mathcal{G},H]\right]+\ldots$

and the decoupling condition

$QH_{\mathrm{eff}}P=0.$ (37)

We partition the original Hamiltonian into an exactly-solvable zero-order part
$H_{0}$ and a perturbation $V$

$H=H_{0}+V$ (38)

and consider an expansion of the generator $\mathcal{G}$ and the interaction
$H_{\mathrm{eff}}$ in powers of $V$

$\mathcal{G}=\mathcal{G}^{[1]}+\mathcal{G}^{[2]}+\mathcal{G}^{[3]}+\ldots\hskip
10.00002pt,\hskip
15.00002ptH_{\mathrm{eff}}=H_{\mathrm{eff}}^{[0]}+H_{\mathrm{eff}}^{[1]}+H_{\mathrm{eff}}^{[2]}+\ldots$
(39)

For convenience, we define the partial sum of the series up to order $n$ as

$\mathscr{G}^{[n]}\equiv\sum_{m=1}^{n}\mathcal{G}^{[m]}.$ (40)

The $n$th order contribution to the effective Hamiltonian is then

$H_{\mathrm{eff}}^{[n]}=\left(e^{\mathscr{G}^{[n]}}He^{-\mathscr{G}^{[n]}}\right)^{[n]}.$
(41)

One can easily verify that $\mathcal{G}^{[n]}$ only contributes to a single
term in Eq. (41). Peeling this term off we have

$H_{\mathrm{eff}}^{[n]}=[\mathcal{G}^{[n]},H_{0}]+\left(e^{\mathscr{G}^{[n-1]}}He^{-\mathscr{G}^{[n-1]}}\right)^{[n]}.$
(42)

Enforcing the decoupling condition (37) yields an equation for
$\mathcal{G}^{[n]}$ in terms of lower-order contributions

$Q[H_{0},\mathcal{G}^{[n]}]P=Q\left(e^{\mathscr{G}^{[n-1]}}He^{-\mathscr{G}^{[n-1]}}\right)^{[n]}P.$
(43)

Equation (43) is of the general form

$[H_{0},X]=Y\,.$ (44)

If we work in the eigenbasis of $H_{0}$ such that
$H_{0}|\phi_{i}\rangle=\epsilon_{i}|\phi_{i}\rangle$, the commutator can be
easily evaluated in terms of the unperturbed energies:

$\langle\phi_{i}|[H_{0},X]|\phi_{j}\rangle=(\epsilon_{i}-\epsilon_{j})\langle\phi_{j}|X|\phi_{j}\rangle\equiv\Delta_{ij}\langle\phi_{j}|X|\phi_{j}\rangle.$
(45)

This suggests that a solution to Eq. (44) can be written as

$X=\frac{Y}{\Delta}+Z\,,$ (46)

where $Z$ is some arbitrary function which commutes with $H_{0}$, and the
superoperator notation introduced in Eq. (19) is used for brevity.

Following this line of reasoning we solve Eq. (43) as777The additional term
like $Z$ in Eq. (46), which commutes with $H_{0}$, will vanish when sandwiched
between $Q$ and $P$, so we need not include that term here.

$Q\mathcal{G}^{[n]}P=Q\frac{\left(e^{\mathscr{G}^{[n-1]}}He^{-\mathscr{G}^{[n-1]}}\right)^{[n]}}{\Delta}P.$
(47)

As we can see, the decoupling condition applies to the $Q\mathcal{G}P$ block
of the generator $\mathcal{G}$, and we have some freedom to choose the rest of
$\mathcal{G}$, namely $P\mathcal{G}P$, $P\mathcal{G}Q$, $Q\mathcal{G}Q$. The
various choices, which we will outline below, result in different effective
Hamiltonians. There are a few important consequences of these choices.

First, in order for the transformation to be unitary, we need the generator to
be anti-Hermitian: $\mathcal{G}=-\mathcal{G}^{\dagger}$. Consequently, the
popular choice $P\mathcal{G}Q=0$ will result in a non-unitary transformation,
and a non-Hermitian effective Hamiltonian. All else being equal, a Hermitian
effective Hamiltonian is preferable, but the significant simplifications that
come with taking $P\mathcal{G}Q=0$ can make this choice attractive.

Figure 4: A schematic illustration showing that the same operator (indicated
by the arrow) can connect _Left:_ a $P$ space configuration to a $Q$ space
configuration as well as _Right:_ a $Q$ space configuration to a $Q$ space
configuration.

Second, the choice $Q\mathcal{G}Q=0$ cannot be enforced in a Fock-space
formulation, and so this choice can only be made when working directly in the
$A$-body Hilbert space. To understand this point 888To our knowledge, this
point was first made in the chemistry literature by Kutzelnigg [180, 181], but
it has not been explicitly stated in the nuclear physics literature. consider
a one-body Fock-space operator which excites a particle from the valence space
to the excluded space. As illustrated in Fig. 4, if this operator acts on a
configuration which belongs to the $P$ space, it will generate a configuration
which belongs to the $Q$ space. However, if that same operator acts on a
$Q$-space configuration which already has some other particle-hole excitation,
it will generate a distinct $Q$-space configuration. So the operator also
connects $Q$-space configurations to $Q$-space configurations. In order to
enforce $Q\mathcal{G}Q=0$ while allowing $Q\mathcal{G}P\neq 0$, the operator
$\mathcal{G}$ needs to ensure that it acts only on $P$-space states, which
means $\mathcal{G}$ must be an $A$-body operator.

Third, if we choose both $P\mathcal{G}P=0$ and $P\mathcal{G}Q=0$, then we have

$Pe^{\mathcal{G}}=P(1+\mathcal{G}+\tfrac{1}{2}\mathcal{G}^{2}+\ldots)=P.$ (48)

If we denote the eigenstate of the full Hamiltonian by $|\Psi_{i}\rangle$ and
the corresponding eigenstate of the effective Hamiltonian $|\psi_{i}\rangle$,
these are related by the similarity transformation

$|\Psi_{i}\rangle=e^{-\mathcal{G}}|\psi_{i}\rangle,$ (49)

and since $|\psi_{i}\rangle$ exists entirely in the $P$ space, Eq. (48)
implies that

$|\psi_{i}\rangle=P|\Psi_{i}\rangle.$ (50)

Physically, the eigenstate of the effective Hamiltonian is just given by the
projection of the full-space wave function to the $P$ space. This means that
$|\Psi_{i}\rangle$ and $|\psi_{i}\rangle$ cannot be simultaneously normalized
to one, and we must employ an “intermediate normalization”
$\langle\Psi_{i}|\psi_{i}\rangle=1.$ On the other hand, if we do not require
$P\mathcal{G}Q=0$, then the above argument no longer holds and the eigenstates
of the effective Hamiltoninan will in general not be simply projections of the
full-space eigenstates.

Finally, before investigating various choices of $\mathcal{G}$, we consider an
iterative method for summing the perturbative series to all orders. To do
this, we notice that selecting out the $n$th order contribution on the right-
hand side of Eq. (47) quickly leads to complicated formulas (see, e.g.,
section III of Ref. [170]). Overall, the right-hand side of Eq. (47) is of
order $n$ and higher, since the lower order terms in $\mathcal{G}$ have been
selected to eliminate the undesired components of $H_{\mathrm{eff}}$ to their
respective orders. This means that if we do not specifically select the $n$th
order terms on the right hand side, but instead take everything, we obtain a
contribution $\mathcal{G}^{[n]}$ which suppresses the $n$th order term in
$H_{\mathrm{eff}}$, as well as some contribution from higher order terms.
These higher order terms will be suppressed during later iterations. Because
this is no longer a strict order-by-order perturbative expansion, we use a
subscript to denote the $n$th iteration, to distinguish it from the
superscript indicating the $n$th-order contribution:

$Q\mathcal{G}_{n}P=Q\frac{e^{\mathscr{G}_{n-1}}He^{-\mathscr{G}_{n-1}}}{\Delta}P.$
(51)

Here, we have defined $\mathscr{G}_{n}\equiv\sum_{m=1}^{n}\mathcal{G}_{m}$ in
analogy with Eq. (40). Defining the transformed Hamiltonian after $n$
iterations as $H_{n}\equiv e^{\mathscr{G}_{n}}He^{-\mathscr{G}_{n}}$, this may
be written as

$Q\mathcal{G}_{n}P=Q\frac{H_{n-1}}{\Delta}P.$ (52)

or

$Q\mathscr{G}_{n}P=Q\mathscr{G}_{n-1}P+Q\frac{H_{n-1}}{\Delta}P.$ (53)

Beginning with $\mathscr{G}_{0}=0$ and iterating Eq. (53) successively
eliminates contributions to $QH_{\mathrm{eff}}P$ of increasing powers of the
perturbation $V$, and for $n\rightarrow\infty$ yields the exact generator:
$\mathscr{G}_{\infty}=\mathcal{G}$.

In the following subsections, we consider choices found in the literature for
fully specifying $\mathcal{G}$, and the consequences of these choices.

####  3.1.2 Lee-Suzuki

We begin with the most restrictive combination
$Q\mathcal{G}Q=P\mathcal{G}P=P\mathcal{G}Q=0$, which allows for the greatest
simplification. Following the notation of Suzuki and Lee [157], for this
choice we write the generator as $\mathcal{G}=-\omega$. The effective
Hamiltonian is

$H_{\mathrm{eff}}=e^{-\omega}He^{\omega}$ (54)

and we have the great simplification that $\omega^{2}=0$ so that
$e^{\omega}=1+\omega$. As we remarked above, this choice will yield a non-
Hermitian effective Hamiltonian which is necessarily formulated in the
$A$-body Hilbert space. The decoupling condition is

$\displaystyle QH_{\mathrm{eff}}P$ $\displaystyle=(1-\omega)H(1+\omega)$ (55)
$\displaystyle=QVP+QH\omega P-Q\omega H_{0}P-Q\omega VP-Q\omega V\omega P=0,$

where in the second line we have used $H=H_{0}+V$. If we take a degenerate $P$
space with energy $\epsilon$ so that $Q\omega PH_{0}P=\epsilon Q\omega P$, we
may rearrage to solve for $\omega$

$\omega=\frac{1}{\epsilon-QHQ}\left(QVP-\omega(PVP+PVQ\omega)\right).$ (56)

As shown by Suzuki and Lee [157], defining

$H_{\mathrm{eff}}=PH_{0}P+PVP+PVQ\omega$ (57)

and iteratively inserting Eq. (56) into itself yields the $Q$-box folded
diagram expansion (7) or, with a different iteration scheme,

$R_{n}=\frac{1}{1-\hat{Q}_{1}-\sum\limits_{m=2}^{n-1}\hat{Q}_{m}\prod\limits_{k=n-m+1}^{n-1}R_{k}}\hat{Q}$
(58)

where $V_{\mathrm{eff}}=R_{\infty}$, $\hat{Q}$ is the same $\hat{Q}$-box
defined in Eq. (6), and
$\hat{Q}_{m}\equiv\frac{d^{m}}{d\epsilon^{m}}(\hat{Q}(\epsilon)$.

As we have pointed out in section 3.1.1, the requirement $Q\omega Q=0$ implies
that this approach must be formulated in the $A$-body space. However, as we
see next, a Fock-space formulation yields the same effective interaction.

####  3.1.3 Generalized Bloch equation

If we drop the requirement $Q\mathcal{G}Q=0$, we have
$P\mathcal{G}Q=P\mathcal{G}P=0$. Following convention, we express the
transformation in terms of the Møller wave operator and its inverse999We use
$\Omega$ for the Møller wave operator for consistency with existing literature
in the present subsection, but caution that it should not be confused with the
Magnus operator used in the IMSRG.

$\Omega=e^{-\mathcal{G}}\hskip 10.00002pt,\hskip
10.00002pt\Omega^{-1}=e^{\mathcal{G}}.$ (59)

Noting that $Q\Omega^{-1}(1-\Omega P)=Q\Omega^{-1}$, we can satisfy the
decoupling condition (37) if

$(1-\Omega P)H\Omega P=0.$ (60)

Again, splitting up $H=H_{0}+V$ and rearranging, this may be written as

$Q[\Omega,H_{0}]P=QV\Omega P-Q\Omega PV\Omega P$ (61)

which is the generalized Bloch equation [184].

The effective interaction in the $P$ space is given by

$PH_{\mathrm{eff}}P=PH\Omega P=PH_{0}P+PV\Omega P.$ (62)

Expanding Eqs. (61) and (62) in powers of $V$ yields a linked expansion for
the Rayleigh-Schrödinger perturbation series [156, 23]. In the nuclear case,
the order-by-order convergence of this series is questionable (cf. Sec. 1).

The second term on the right hand side of Eq. (61) can be represented by
“folded” diagrams [156, 23]. The $\hat{Q}$-box method described in section 2.1
amounts to a perturbative expansion of the first term in Eq. (61), followed by
a summation of certain higher-order terms in the folded diagram series. We
note that if we express the wave operator as $\Omega=P+\chi$, where $\chi$ is
called the “correlation operator”, then Eq. (61) is equivalent to Eq. (55)
with $\chi=\omega$. Interestingly, while the Lee-Suzuki approach required
$Q\omega Q=0$, and therefore could not be expressed in a Fock-space
formulation, the Bloch equation approach does not use that constraint and so
may be formulated in Fock-space. Evidently, the Bloch equation approach does
not make any reference to $Q\mathcal{G}Q$, and so setting it to zero does not
alter the resulting effective Hamiltonian. For a discussion of the differences
between the Lee-Suzuki scheme and the $\hat{Q}$-box approach, see e.g. Ref.
[28].

####  3.1.4 Shell model coupled cluster

Instead of working with the wave operator, we may directly work with the
generator $\mathcal{G}$. Following the notation of Ref. [179], we write
$\mathcal{G}$ as $-S$. We employ the iterative procedure laid out in section
3.1.1 to obtain an iterative expression for $S$:

$S_{n}=S_{n-1}-Q\frac{H_{n-1}}{\Delta}P$ (63)

where we have denoted the similarity-transformed Hamiltonian101010In practical
applications, the procedure is applied to the CC Hamiltonian (33), (c.f. Sec.
2.4) at the $n$th step $H_{n}=e^{-S_{n}}He^{S_{n}}$. The effective Hamiltonian
is then $H_{\mathrm{eff}}=H_{\infty}$. In arriving at Eq. (63), we have
implicitly assumed that the terms in $QSQ$ are only those which also
contribute to $QSP$. Consequently, if a term in $S$ only connects $Q$
configurations, then it is taken to be zero.

Often, the iterations are better behaved with the help of a convergence
factor, which we denote $dt$ to connect with the formulation of Ref. [179].
Multiplying the second term of Eq. (63) and taking the limit $dt\rightarrow
0$, we can reinterpret it as a flow equation:

$\frac{dS}{dt}=-\eta(t)\equiv-Q\frac{H(t)}{\Delta}P.$ (64)

Taking the $P$ space to be a single Slater determinant, Eq. (63) yields an
iteration scheme to solve the coupled cluster equations, while defining the
$P$ space in terms of a valence space yields the shell model coupled cluster
approach [179] described in section 3.1.4. As discussed above, the requirement
$PSQ=0$ means that $H_{\mathrm{eff}}$ is not Hermitian in this approach.
Additionally, we see that the SMCC effective interaction is equivalent to the
other two non-Hermitian effective interactions discussed in sections 3.1.2 and
3.1.3, as long as no approximations are made.

####  3.1.5 Canonical perturbation theory

We next drop the condition $Q\mathcal{G}P=0$, enabling us to enforce
$\mathcal{G}^{\dagger}=-\mathcal{G}$ so that the transformation is unitary and
$H_{\mathrm{eff}}$ is Hermitian. We retain the condition $Q\mathcal{G}Q=0$,
with the consequence that we cannot express the theory in terms of Fock-space
operators. However, this restriction greatly simplifies the analysis. For
consistency with the literature, we write the generator as $\mathcal{G}=-G$.
It can be shown [178, 185] that the operator $G$ is related to the operator
$\omega$ from section 3.1.2 by111111Here, the hyperbolic arctangent of an
operator is defined in terms of its Taylor series expansion [178].

$G=\mathrm{arctanh}\left(\omega-\omega^{\dagger}\right)$ (65)

and the transformation is [185]

$e^{G}=(1+\omega-\omega^{\dagger})(1+\omega\omega^{\dagger}+\omega^{\dagger}\omega)^{-1/2}.$
(66)

The effective Hamiltonian obtained in this approach is the Hermitized version
of the effective Hamiltonian resulting from Eq. (55). To connect this result
with the OLS approach presented in section 2.2, we write that transformation
out explicitly:

$H_{\mathrm{eff}}=\frac{\mathcal{U}_{P}^{\dagger}}{\sqrt{\mathcal{U}_{P}^{\dagger}\mathcal{U}_{P}}}\mathcal{U}H\mathcal{U}^{\dagger}\frac{\mathcal{U}_{P}}{\sqrt{\mathcal{U}_{P}^{\dagger}\mathcal{U}_{P}}}.$
(67)

By inserting a sum over the eigenstates of $H$, and using
$|\Psi_{i}\rangle=(1+\omega)P|\Psi_{i}\rangle$, one can show [162, 163] that
$\mathcal{U}^{\dagger}\mathcal{U}_{P}=(1+\omega)P$, and the OLS effective
interaction is equivalent to $Pe^{-G}He^{G}P$.

An iterative scheme very similar to the one described in section 3.1.1 was
proposed by Suzuki [171], but not directly pursued further. The unitary model
operator approach (UMOA) [185] follows this formalism, with the valence
cluster expansion (see section 3.2.2) carried out on the generator, rather
than on the effective Hamiltonian. So far, studies with the UMOA have focused
on ground state energies of closed-shell nuclei [186], so we will not discuss
it further here.

####  3.1.6 Fock-space canonical perturbation theory

If we desire a Hermitian effective operator with a Fock-space decomposition,
then we should drop the restriction $Q\mathcal{G}Q=0$, leaving only
$P\mathcal{G}P=0$. Writing out Eq. (47) order-by-order with the requirement
$\mathcal{G}^{\dagger}=-\mathcal{G}$ yields the canonical perturbation theory
of Primas [187] and Klein [170]. Interestingly [181], the resulting expansion
is different from the expansion obtained in section 3.1.5, i.e. the Fock-space
and $A$-body space formulations are not equivalent, in contrast to what was
found for the non-Hermitian formulation. This approach has not been pursued in
the nuclear physics literature.

####  3.1.7 Unitary coupled cluster

Alternatively, we may take $P\mathcal{G}P=0$ and
$\mathcal{G}^{\dagger}=-\mathcal{G}$ and follow the iterative procedure of
section 3.1.1 to obtain

$\mathcal{G}_{n}=\mathcal{G}_{n-1}+Q\frac{H_{n-1}}{\Delta}P+P\frac{H_{n-1}}{\Delta}Q.$
(68)

This yields a unitary coupled cluster expansion for the effective interaction.
As with the SMCC solution, this may be recast as a differential equation (here
we use $s$ instead of $t$)

$\frac{d\mathcal{G}(s)}{ds}=Q\frac{H(s)}{\Delta}P+P\frac{H(s)}{\Delta}Q.$ (69)

This approach has also not been persued in nuclear physics, although it is
very closely related to the Magnus formulation of the IMSRG, as we will show
below.

####  3.1.8 Magnus IMSRG

Finally, we may drop the constraint $P\mathcal{G}P=0$ and instead specify
$\mathcal{G}$ by the requirement that it should reproduce the flowing
Hamiltonian $H(s)$ along its entire trajectory. Following Ref. [174], we write
$\mathcal{G}=\Omega(s)$, which we call the Magnus operator, and require

$e^{\Omega(s)}He^{-\Omega(s)}=H(s)$ (70)

where $H(s)$ is the solution of the flow equation (11). This may be
reorganized as a flow equation for the Magnus operator $\Omega(s)$, as
described in section 2.3.2. Considering the first few terms in the series, we
have

$\frac{d\Omega(s)}{ds}=\eta(s)-\frac{1}{2}[\Omega(s),\eta(s)]+\ldots$ (71)

If we choose the White generator (17), which may be written as

$\eta^{\mathrm{Wh}}(s)=Q\frac{H(s)}{\Delta}P+P\frac{H(s)}{\Delta}Q,$ (72)

and neglect all terms on the right hand side of Eq. (71) aside from the first
one, we recover the unitary coupled cluster equation (69). The difference
between Magnus IMSRG and unitary coupled cluster then lies in the commutator
terms of Eq. (71). A perturbative analysis reveals that the leading-order
effect of the first commutator term is to induce contributions to $P\Omega P$
and $Q\Omega Q$ at third order. If the transformation is evaluated exactly,
these terms of course have no effect on the resulting observables. However, if
approximations are made—as they inevitably must be—then these terms may
produce a different result. This has not yet been investigated in detail.

Integrating Eq. (71) numerically with a step size $ds=1$ —again neglecting all
but the first term— we find, following the discussion leading to Eq. (53),
that the first integration step yields a generator which satisfies the
decoupling condition to first order in perturbation theory (cf. Appendix C).
Likewise, the second step in $ds$ satisfies the decoupling condition to second
order, and the $n$th step satisfies decoupling to $n$th order. Thus, numerical
integration of the flow equation with step size $ds=1$ corresponds to an
order-by-order summation of the perturbation series. If instead we take a
smaller step size, $ds=0.5$, then after the first integration step we will
have only suppressed half of the first-order term in the decoupling condition.
After the second integration step, taking us to $s=1$, we have suppressed half
of the _remaining_ first order term, as well as half of the second order term.
Taking the continuous limit $ds\rightarrow 0$, we find the first-order off-
diagonal piece suppressed as $e^{-s}$, with the higher-order terms also
suppressed at the same rate.

In light of this discussion, we can view the numerical integration of the flow
equation (71) with some finite step-size $ds$ as a summation of the
perturbative expansion [145, 146], with the step size specifying anything from
an order-by-order summation ($ds=1)$ to all orders at once ($ds\rightarrow
0$). This connection may have important consequences in cases where the
perturbative expansion does not converge (see section 6.2).

###  3.2 Approximation schemes

An exact evaluation of the formulas for $H_{\mathrm{eff}}$ presented in the
preceding subsections will inevitably be at least as expensive as a direct
diagonalization of the Hamiltonian in the full Hilbert space—precisely the
task we set out to avoid. The utility of the effective interaction framework
is that it facilitates approximations which greatly reduce the required effort
while minimally impacting the accuracy of the computed quantities of interest,
namely observables related to low-lying eigenstates. Within the shell-model
context, this implies the need for some sort of cluster truncation.

Perhaps the most straightforward approximation scheme is a truncation in in
powers of the residual interaction $V$, i.e. perturbation theory.
Unfortunately, in nuclear physics the effective interaction often converges
slowly in powers of $V$, and as discussed in section 6.2, the intruder-state
problem suggests that in most cases the perturbation series is divergent.

Nonperturbative truncation schemes have been made essentially along two lines:
either a cluster truncation is imposed within a Fock-space formulation, as for
IMSRG and CC (cf. Secs. 2.3.2 and 3.1.4), or the problem is solved directly in
the $A$-body system for a few valence particles, followed by a cluster
expansion.

####  3.2.1 Fock-space cluster truncation

When working in a Fock-space formulation, it is natural to perform a cluster
truncation on the generator $\mathcal{G}$ or on all operators, typically
limiting them to consist of zero-, one- and two-body pieces. From a practical
point of view, such a truncation is a necessity; keeping three-body terms is
unpleasant but feasible, while the need for e.g. six-body terms would be
sufficiently onerous to render the method useless.

For the special case of a one-dimensional $P$ space (i.e. a single-reference
calculation), using the non-hermitian formulation of section 3.1.4, truncating
$\mathcal{G}$ to one- and two-body operators is equivalent to coupled cluster
with singles and doubles (CCSD) [25]. Here, we find the desirable feature that
the Baker-Campbell-Hausdorff expansion (36) formally truncates after a finite
number of nested commutators (four in the CCSD approximation, if $H$ has at
most two-body terms).

For the case of interest in the context of the shell model, with the dimension
of the $P$ space greater than one, the Baker-Campbell-Hausdorff expansion does
not formally truncate [179]. One approach to this issue is to truncate the
series at a finite order of perturbation theory, or else at a finite power of
$\mathcal{G}$ (see, e.g. [177]). Alternatively, one can specify a form for the
Fock-space operators, e.g. retaining one- and two-body terms while discarding
the rest, allowing the series to be evaluated iteratively [169, 175]. While
the series remains infinite with this truncation, in most cases of interest it
is found that the series converges so that for a given precision only a finite
number of nested commutators must be evaluated. Importantly, this truncation
scheme retains only connected diagrams, and so maintains size extensivity.
This approach is used in the VS-IMSRG and SMCC described above, and in the
canonical transformation theory of Yanai et al [176, 188].

When operators are normal ordered with respect to a finite-density reference,
many-body operators can feed back into fewer-body operators through the
commutators in the Baker-Campbell-Hausdorff expansion. However, the reduction
in particle-rank of an operator always comes with an occupation number (see
the flow equation in Appendix A), corresponding to a factor of the density,
and so the discussion in Sec. 2 about the cluster hierarchy justifies this
truncation.

####  3.2.2 Valence cluster expansion

The other approximation scheme is to work within the $A$-body Hilbert space
formulation and build up the effective interaction in order of increasing
cluster rank [163]. One diagonalizes the
$A_{\mathrm{core}},A_{\mathrm{core}}+1,A_{\mathrm{core}}+2$ systems
successively, and extracts the consistent core energy, single-particle
energies, and two-body matrix elements by subtracting the contributions of
lower particle rank. One could continue this procedure to obtain higher-body
effective interaction, with rapidly increasing effort. Instead, assuming that
the effective interaction has a sufficiently convergent cluster expansion, the
effective interaction obtained in the zero-, one-, and two-valence-particle
systems can then be applied to systems with more valence particles. For self-
bound systems like nuclei, one must take care to properly treat the mass-
dependence of the intrinsic kinetic energy in the construction procedure for
the effective interaction [163, 165, 166], although this effect becomes less
important for heavier systems.

This scheme is used in the Okubo-Lee-Suzuki approaches based on coupled
cluster (CCEI) [140, 166] and the no-core shell model [162, 163, 165], as well
as the $\hat{Q}$-box approach of Kuo and collaborators [26, 158, 30]. One
potential drawback of the valence cluster expansion is that the optimal
cluster decomposition for two valence particles might be significantly
different from the optimal decomposition for many valence particles.
Considering the $sd$ shell as a specific example, the appropriate mean field
for an empty valence space, corresponding to 16O, will be quite different from
the appropriate mean field for a filled valence space, corresponding to 40Ca,
and so one would expect that different single-particle energies would be
optimal.

Another potentially more serious issue arises from a perturbative analysis. As
was discussed in Sec. 3.1.1, the requirement $Q\mathcal{G}Q=0$ cannot be
enforced in a Fock-space formulation. If it is enforced in the $A$-body
formulation, an analysis reveals that disconnected diagrams arise [189, 181].
For example, in a system with four valence particles and a purely two-body
interaction, a disconnected four-body contribution to the effective
interaction arises at fourth order. Such a contribution — essentially two-body
interactions between two independent pairs — is not subject to the arguments
at the beginning of Sec. 2 about short-range interactions at low density,
because it does not depend on the separation between the clusters. We may
expect such terms to be combinatorially enhanced, which would be a serious
problem — this calls for further investigation.

The above point may provide some explanation for the surprising finding in
Ref. [166], where a single interaction for the $sd$ shell obtained with the
CCEI method produced impressive agreement with experimental binding energies
throughout the shell. While at first glance such a result is cause for
celebration, a closer look suggests trouble. As a specific example, other _ab
initio_ calculations [190, 126, 191, 143] (including coupled cluster) using
the same input interaction find 40Ca to be over-bound by nearly 40 MeV, with
relatively small variation among the calculations, while the CCEI result is
under-bound by just 2 MeV. The supplemental material of Ref. [166] makes note
of this, since it is unreasonable to expect that the CCEI method should be
more accurate than the coupled cluster method upon which it is based. In the
present context, we may speculate that the combination of _(i)_ the over-
binding inherent in the input force, _(ii)_ missing valence many-body effects
(see Sec. 4.2) enhanced by the effect of disconnected diagrams, and _(iii)_
truncation errors in high-lying eigenvalues from the equation-of-motion
coupled cluster method [192] incidentally conspire to cancel out in the $sd$
shell. However, one can and should not rely on such a cancellation in general.

##  4 Three-body forces and the connection with phenomenological adjustments

Shortly after Yukawa’s formulation of the nuclear interaction in terms of pion
exchange, it was pointed out [193] that a description of a quantum field
theory in terms of an instantaneous (or, equivalently, energy-independent)
potential inevitably leads to three-body and higher-body forces. The
connection between these many-body forces and nuclear saturation was also
suspected early on [194, 195, 70], although the calculations were necessarily
schematic. Likewise, mean-field calculations using a Skyrme [196, 197] or
Gogny [198] parameterization of the force require a three-body, or density-
dependent two-body, term.

Figure 5: Three-body forces generated from _Left:_ elimination of the $\Delta$
isobar degree of freedom, and _Right:_ elmination of an excitation to a
$Q$-space configuration.

Of course, even if the initial interaction were solely of a two-body nature,
the effective interaction in the valence space will still in general contain
three-body and higher-body forces. In fact, these “effective” or “induced” 3N
forces and the “genuine” 3N forces are essentially of the same origin: the
elimination of degrees of freedom. This is illustrated in Fig. 5, where the
elmination of the $\Delta$ isobar degree of freedom and the elimination of an
excitation toa $Q$-space configuration both lead to effective 3N interactions.
Some previous studies have found the induced and “genuine” 3N forces are of
comparable magnitude [199, 200], although this is depends on the
renormalization scheme and scale of the interaction.

As a practical matter—due to the ambiguity in producing a three-body term
consistent with the two-body interaction, as well as the difficulty in
handling a three-body term in a many-body calculation—explicit three-body
forces have historically been neglected in shell model treatments, although
there were some exploratory calculations (see e.g. [201, 202, 203, 199, 200]).
There have also been more recent calculations evaluating 3N forces in the
valence space either perturbatively [204] or explicitly [179]. Zuker and
collaborators [205, 18] argued that the main effect of the three-body forces
should be to modify the monopole (i.e. diagonal, $J$-averaged) component of
the effective interaction. This argument simultaneously justified the omission
of explicit three-body terms as well as the phenomenological adjustment of
monopole terms in the effective interaction, which resulted in excellent
reproduction of the experimental data. Further supporting this point of view
were _(i)_ the fact that the various realistic NN interactions produced
similar shell model matrix elements121212This observation is easily understood
from the RG/EFT point of view—the various potentials differ in their high-
momentum content but reproduce the same low-momentum physics [206]., allowing
little room for improvement, and _(ii)_ the observed improvement in
spectroscopy of light nuclei obained of quantum Monte Carlo and no-core shell
model calculations when explicit three-body forces were included [207, 208].

An important demonstration of the effect of three-body forces in the shell
model was a calculation showing that three-body forces could help explain the
location of the neutron dripline in oxygen [132], followed by an explanation
of the $N=28$ magic number in the calcium isotopes [134]. These calculations
used a normal-ordering approximation (see Sec. 2.3.1) for the three-body force
and obtained essentially the monopole effect described by Zuker, although they
used empirical single-particle energies and scaled the two-body matrix
elements by $A^{1/3}$ as in phenomonological calculations. The same effect was
soon confirmed in _ab initio_ calculations without phenomenological
adjustments [209, 125, 210].

The first VS-IMSRG calculations of the oxygen isotopes did not obtain the
correct dripline [139], even though three-body forces were included in the
normal-ordered approximation. More troubling, the heavier oxygen isotopes were
systematically overbound by approximately 10 MeV. The issue was that the
normal ordering of the Hamiltonian used the core wave function as a reference
state in these initial VS-IMSRG calculations, hence the effects of three-body
interactions between valence nucleons were not properly captured. This
deficiency was remedied by the use of ensemble normal ordering (ENO) [143],
which enables an approximate treatment of the effect of three-body forces that
does not degrade as valence particles are added. This echoes the results of
previous investigations of the effects of three-body forces in the shell model
[211, 200, 199, 29]. Since only a brief account of ENO has been given in the
literature [143], we provide a more detailed description in the following
section.

###  4.1 Ensemble normal ordering

When performing a VS-IMSRG calculation, a natural choice for the normal
ordering reference $|\Phi\rangle$ is the core of the valence space. This
allows an approximate treatment of 3N forces in which a sum over particles in
the core yields effective one-body and two-body forces in the valence space.
However, VS-IMSRG calculations of the oxygen isotopic chain using chiral NN+3N
forces overpredicted the binding energy of neutron-rich oxygen nuclei compared
to an earlier MR-IMSRG study with the same interactions [139, 125].
Calculations involving both protons and neutrons in the valence space yielded
even more significant overbinding [142].

This discrepancy was essentially due to the fact that the normal ordering in
the MR-IMSRG calculation is performed directly with respect to the system of
interest, not with respect to the core of the valence space. This meant that,
in 24O for example, the MR-IMSRG was better capturing the 3N interactions
between the 8 valence neutrons. Indeed, taking the normal ordering reference
to be the nearest closed-shell nucleus brought the VS-IMSRG binding energies
back in line with the MR-IMSRG results [142].

Figure 6: (a) The binding energy as a function of mass number $A$ for the
oxygen isotopes, calculated with a variety of many-body methods[125, 212, 140,
126, 213]. The blue curve labeled IMSRG(SM) corresponds to the scheme in (b)
with the core taken as the normal-ordering reference, while the red curve
labeled IMSRG(ENO) corresponds to (c) with an ensemble reference. Adapted from
[143].

This approach was then generalized to treat systems that are not close to any
closed sub-shell, by allowing fractional occupation numbers. As an example,
consider 19O, which in a naive shell model picture has three neutrons in the
$0d_{5/2}$ orbit on top of a closed 16O core. Equivalently, it could be
considered as three neutron holes in 22O. When using an 16O reference, the
occupation number for the neutron $0d_{5/2}$ orbit would be 0, while it would
be 1 for an 22O reference. The 16O reference will underestimate the missing
three-body effects, while the 22O reference will overestimate them. The
compromise is then to take the occupation to be 0.5, i.e. filling the orbit
half way. This strategy is frequently used in mean-field theory, and it is
known as the “equal-filling approximation” or simply the “filling
approximation”. (For an application in chemistry, see [188]).

The question then arises: what reference state (if any) is actually being used
when we select fractional occupation numbers? As explained in Ref. [214, 215]
the equal filling approximation can be framed in terms of a mixed state or
ensemble, in the sense of quantum statistical mechanics, specified by a
density matrix131313Not to be confused with the one-body density matrix
$\rho_{pq}=\langle\Psi|a^{\dagger}_{p}a_{q}|\Psi\rangle$.

$\rho=\sum_{\alpha}c_{\alpha}|\Phi_{\alpha}\rangle\langle\Phi_{\alpha}|$ (73)

for some set of coefficients $c_{\alpha}$. (Here $\alpha$ labels different
Slater determinants). The expectation value of an operator $\mathcal{O}$ in
the ensemble is obtained by a trace over the density matrix:
$\langle\mathcal{O}\rangle=\mathrm{Tr}\left[\mathcal{O}\rho\right]=\sum_{\alpha}c_{\alpha}\langle\Phi_{\alpha}|\mathcal{O}|\Phi_{\alpha}\rangle$.

As discussed in Sec. 2.3, for a single reference $|\Phi\rangle$, the normal
order of a pair of creation and annihilation operators is the one which gives
zero expectation value in the reference. The normal order of a string of more
than two creation/annihilation operators can be chosen so that every pair of
operators in the string is in normal order.

Wick’s theorem [172], which expresses a string of creation/annihilation
operators in terms of normal-ordered strings and contractions, was extended to
more general reference states by Kutzelnigg and Mukherjee [173]. In this case,
the normal ordering is still defined so that the reference expectation value
of a normal ordered string of creation/annihilation operators vanishes.
However, the concept of a contraction now involves one-body, two-body, and
higher-body density matrices, which encode the correlations contained in the
reference. This generalized Wick’s theorem is used, for instance, to formulate
the MR-IMSRG [125, 146].

As shown by Kutzelnigg and Mukherjee, Wick’s theorem also applies to a mixed-
state, or ensemble, reference:

$\langle\\{\\!a^{\dagger}_{p}a_{q}\\}\rangle=\mathrm{Tr}\left[\\{\\!a^{\dagger}_{p}a_{q}\\}\rho\right]=\sum_{\alpha}c_{\alpha}\langle\Phi_{\alpha}|\\{\\!a^{\dagger}_{p}a_{q}\\}|\Phi_{\alpha}\rangle=0.$
(74)

This is an extension of the original formulation of the finite-temperature
Wick’s theorem [216, 217, 218], which only applied to the _expectation value_
of an operator in the ensemble. The formulation of Kutzelnigg and Mukherjee,
on the other hand, is an operator identity, just like the original zero-
temperature Wick’s theorem.

Our goal is then to find an ensemble such that contractions have the same form
as in the single-reference case, except that the occupation of an orbit may
have some non-integer value between 0 and 1. That is, we want all two-body and
higher-body irreducible density matrices, as well as the off-diagonal one-body
density matrix, to vanish. Such an ensemble must necessarily contain a
variable number of particles. This may be easily understood by considering a
single particle placed in two levels. If we require that there always be a
fixed total number of particles (like in canonical or micro-canonical
ensembles, for instance), then the occupation of one level implies that the
other level must be empty and so the occupations are correlated, leading to an
irreducible two-body density matrix. A grand canonical ensemble, on the other
hand will meet our needs.

The ensemble that has been used in VS-IMSRG calculations published thus far
corresponds to the zero-temperature limit of a finite-temperature Hartree-Fock
calculation [218], with the chemical potential chosen to fix the average
particle number. To illustrate the application of Wick’s theorem, we consider
a single level with degeneracy $d=2j+1$, so a configuration has $\mathcal{N}$
particles with $0\leq\mathcal{N}\leq 2j+1$. The contraction of two operators
is given by

$\displaystyle\contraction{}{a}{{}^{\dagger}_{p}}{a}a^{\dagger}_{p}a_{q}=\langle
a^{\dagger}_{p}a_{q}\rangle$
$\displaystyle=\delta_{pq}\frac{1}{\mathcal{Z}}\sum_{\mathcal{N}=0}^{2j+1}\begin{pmatrix}2j\\\
\mathcal{N}-1\end{pmatrix}e^{\beta(\epsilon-\mu)\mathcal{N}}$ (75)
$\displaystyle=\delta_{pq}\frac{1}{\mathcal{Z}}\sum_{\mathcal{N}=0}^{2j+1}\begin{pmatrix}2j+1\\\
\mathcal{N}\end{pmatrix}\frac{\mathcal{N}}{2j+1}e^{\beta(\epsilon-\mu)\mathcal{N}}$
$\displaystyle=\delta_{pq}\frac{1}{2j+1}\langle\mathcal{N}\rangle$
$\displaystyle=\delta_{pq}n_{p}.$

Here $\beta$ is the inverse temperature, $\mu$ is the chemical potential,
$\epsilon$ is the energy of the level, $\mathcal{Z}$ is the partition
function, and the binomial coefficients count how many of the configurations
with $\mathcal{N}$ particles will have orbit $p$ occupied. Of course, other
ensembles may be selected, and they need not be thermal ensembles, i.e.
multiple levels could be fractionally occupied and there need not be a
connection between the energy of a level and its occupation [215]. Another
very reasonable choice of occupations is to use natural orbitals, or a
perturbative approximation of them [147].

In actual calculations, this ensemble need not be explicitly constructed; we
only use the corresponding definition of the contraction when we use Wick’s
theorem. In fact, there is another reference which can be constructed to
produce the same fractional occupations. Instead of employing an ensemble
state, we may use a single-determinant reference built from a single-particle
basis which is slightly different from the one used in the calculation. To
fractionally fill an orbit $p$, we admix in some other “inert” orbit
$\mathcal{Q}$ which is orthogonal to all of the single-particle states used in
our calculation

$\begin{pmatrix}|p\rangle\\\
|\mathcal{Q}\rangle\end{pmatrix}\rightarrow\begin{pmatrix}|\bar{p}\rangle\\\
|\bar{\mathcal{Q}}\rangle\end{pmatrix}=\begin{pmatrix}\sqrt{n}&\sqrt{1-n}\\\
-\sqrt{1-n}&\sqrt{n}\end{pmatrix}\begin{pmatrix}|p\rangle\\\
|\mathcal{Q}\rangle\end{pmatrix}$ (76)

where $0\leq n\leq 1$. If we choose a reference in which the orbit $\bar{p}$
is filled, i.e. $\langle\Phi|a^{\dagger}_{\bar{p}}a_{\bar{p}}|\Phi\rangle=1$,
then the occupation in terms of the original orbit $p$ is
$\langle\Phi|a^{\dagger}_{p}a_{p}|\Phi\rangle=n$. Because the reference $\Phi$
is a single Slater determinant, all higher-body density matrices vanish
automatically. In addition, there will be a non-zero occupation of the inert
orbit $\mathcal{Q}$:
$\langle\Phi|a^{\dagger}_{\mathcal{Q}}a_{\mathcal{Q}}|\Phi\rangle=(1-n)$, as
well as off-diagonal one-body densities
$\langle\Phi|a^{\dagger}_{p}a_{\mathcal{Q}}|\Phi\rangle=\sqrt{n(1-n)}$, which
are not desirable. However, we have asserted that the orbit $\mathcal{Q}$ is
inert. By this we mean that $a^{\dagger}_{\mathcal{Q}}$ and $a_{\mathcal{Q}}$
do not appear in any operator we consider, and we may neglect terms involving
orbit $\mathcal{Q}$ without changing the physics. While the introduction of
inert orbits might seem contrived, it is no more contrived than the ensemble
with variable particle number. Indeed, we could say that the inert orbits live
somewhere in the reservoir that supplies the additional particles.

Again, for practical purposes, it is irrelevant whether the reference is an
ensemble or is constructed with an inert orbit mixed in. What matters is that
we may use Wick’s theorem with fractional occupation numbers, and that this
procedure constitutes an exact rewriting of our operators—given an operator
which is normal ordered with fractional occupations, we can reconstruct the
operator normal-ordered with respect to the true vacuum. The importance of
this point is that by employing fractionally-filled orbitals we have not
introduced an additional approximation. If we retain all the induced operators
up to $A$-body operators, then the IMSRG calculation is exact. What the
fractional filling does is reduce the impact discarding the residual three-
body terms has on the low-lying states. Because standard shell model codes
typically work with valence particles (not valence holes), after the IMSRG
decoupling, we again use Wick’s theorem to rewrite all operators in normal
order with respect to the core (which is a single Slater
determinant)141414This “re-normal-ordering” is easily achieved by using (24),
replacing $n_{a}\rightarrow(n_{a}^{\mathrm{new}}-n_{a}^{\mathrm{old}})$..

Certainly, an uncorrelated ensemble reference is a crude approximation of the
exact wave function, and one might envision that a correlated reference state
as used in the MR-IMSRG could do better [125, 146]. However, as argued in Ref.
[188], it is not clear that this is the best way to proceed in a valence space
context. If correlation effects are included in the reference in order to
better approximate a particular state, then this might well deteriorate the
description of other low-lying states (which contain different correlations),
leading to a worse overall description of the spectroscopy.

###  4.2 Mass dependence of the effective interaction

A significant consequence of the ensemble normal ordering (ENO) procedure is
that a different valence space interaction is obtained for each nucleus. It is
important to emphasize here that because the procedure does not involve any
fitting to data, there is _no loss of predictive power_. The ENO should be
considered as a technique for reducing the impact of the truncation to two-
body operators. In terms of computational effort, the need to generate a new
interaction for each nucleus makes a study of the full $sd$ shell more
laborious, but still manageable. For nuclei in middle of the $pf$ shell, the
exponential scaling of the valence space diagonalization catches up with the
polynomial scaling of the VS-IMSRG and so generating the effective interaction
takes about as long as the shell model calculation that uses it.

The need for some mass-dependence of the effective interaction has been known
for a long time. The $sd$ shell interactions of Kuo and collaborators [6, 7]
yielded a good description of spectroscopy for a few valence particles or
valence holes, but agreement deteriorated for mid-shell systems [219].
Investigations by Chung and Wildenthal [52, 10] suggested that a single
phenomenological adjustment could not remedy the situation, and a scaling of
the two-body matrix elements according to $A^{0.3}$ was introduced. This
prescription has been adopted in many later treatments [11, 134, 30]. The
scaling is typically justified in terms of the increasing nuclear radius
changing the optimal harmonic oscillator frequency [10, 56, 11, 144]. While
such an argument would suggest that the core and single-particle energies
should also change with mass, these effects could in principle be absorbed
into the scaling of the two-body matrix elements [10]. On the other hand, the
need for mass-dependence of two-body matrix elements could be interpreted as a
signal of non-negligible three-body terms in the effective interaction, and
indeed this has been suggested a number of times [211, 199, 200].

We may expect that ensemble normal ordering should capture both the effects of
a changing mean field and of the residual three-body effective
interaction151515Indeed, these effects are not entirely distinct; the induced
three-body interaction depends on the choice of reference. . Fig. 7 displays
binding energies per nucleon obtained for oxygen isotopes and $N=Z$ nuclei in
the $sd$-shell nuclei with the USDB interaction, both with and without the
mass scaling of the TBMEs. These are compared to the binding energies of the
same nuclei calculated using the VS-IMSRG with and without ensemble normal
ordering. It is evident that the ensemble normal ordering has qualitatively
the same effect as the scaling of the TBMEs, athough an investigation of the
VS-IMSRG TBMEs reveals no such smooth scaling (the effect is largely captured
in the core and single-particle energies).

Figure 7: Energy per nucleon for _Left_ the oxygen isotopes $16\leq A\leq 28$,
and _Right_ the $N=Z$ nuclei in the $sd$ shell, obtained with VS-IMSRG using
the EM1.8/2.0 interaction compared with the results obtained with the USDB
interaction. The thinner lines indicate the effect of turning off the ensemble
normal ordering (ENO) in the VS-IMSRG calculation, or turning off the
$A^{0.3}$ scaling of two-body matrix elements in the USDB interaction. Figure
8: (a) Single particle energies in the $sd$ shell obtained with the VS-IMSRG
normal ordered with respect to a 28Si reference, using an NN interaction, with
or without the 3N piece. (b)Neutron-neutron and proton-neutron monopoles of
the two-body matrix elements, with and without the 3N force. (c) The
difference between matrix elements obtained with NN only and NN+3N, with and
without a monopole correction. All calculations use the EM 1.8/2.0 interaction
of Ref. [220].

Figure 8 shows the single-particle energies and monopoles of two-body matrix
elements (TBMEs) obtained for a 28Si reference with and without explicit 3N
forces. Including the 3N force has a significant impact on the single-particle
energies; indeed, the neutron SPEs are shifted closer to the USD [57] values
(USD does not include the Coulomb interaction). In panel (b) of Fig. 8, we see
that the effect of the 3N interaction on the TBME monopoles is repulsive, as
expected from binding energy calculations, and that they shift the monopoles
towards the USDB values. In panel (c) of Fig. 8 we show the difference between
each of the two-body matrix elements obtained with the NN only and NN+3N
interactions (black crosses) and the difference when the NN only monopoles
have been shifted to the NN+3N values. The monopole shift does not yield
perfect agreement—there is still some scatter in the red markers—but the
remaining discrepancy is approximately Gaussian and centered on zero. It is
not unreasonable that there would be moderate cancellation between the
remaining terms, and that the monopole correction would approximately account
for the missing 3N forces, as claimed by Zuker et al. [205]. For a related
approach using density functionals to inform the monopole correction, see Ref.
[221].

Based on our discussion, we can conclude that the “standard” phenomenological
adjustments made to shell model interactions can be understood essentially in
terms of the effect of missing (normal-ordered) three-body forces (see e.g.
the discussion in section 8 of Ref. [156], as well as Refs. [222, 211, 199,
200]).

  * •

Historically, the core energy was taken from experiment (typically no effort
was made to calculate it consistently from the input force), and modern _ab
initio_ calculations have confirmed the importance of three-body forces to
binding energies [223, 209, 190].

  * •

Likewise, single-particle energies were typically taken from experiment, as
the ones obtained from the NN interaction did not reproduce the spectra of
one-valence-particle systems. The normal-ordered contribution of three-body
forces to the single-particle energies essentially accounts for this
discrepancy.

  * •

Even with the core and single-particle energies taken from experiment,
realistic NN forces typically did not give good spectroscopy, and needed
phenomenological adjustment. Zuker [205] argued that the most important
adjustment was of the two-body monopoles, and that this shift should be
understood in terms of missing three-body forces. Indeed, when the three-body
contribution to the normal-ordered two-body interaction is taken into account,
no phenomenological shifts are needed. Moreover, as shown in Figure 8, the
bulk of the discrepancy between an interaction derived from only NN forces and
one including 3N effects can be corrected by a shift of the monopoles. The
remaining discrepancy (the “multipole” terms) is approximately Gaussian and
centered on zero, so that the net effect will be in general small.

  * •

Finally, the $\sim A^{1/3}$ scaling of two-body matrix elements can be
understood as a way to capture the bulk effects of three-body forces among
valence particles. This same physics is captured by employing ensemble normal
ordering (see also [204, 224]).

##  5 Applications

As discussed previously, shell model calculations with _ab initio_
interactions allow us to confront our starting point, the underlying chiral
two- plus three-nucleon force of our choice, with a wealth of available
experimental data. Until recently, applications have been primarily focused on
ground- and excited state energies, with very encouraging results. In the
following, we will highlight selected examples.

###  5.1 Ground and Excited States of $sd-$Shell Nuclei

In Sec. 4, we discussed the importance of chiral three-nucleon forces for the
correct description of nuclear shell structure. Their effect on the location
of the oxygen dripline was one of the first high-profile applications of _ab
initio_ interactions in the nuclear shell model [132], which has subsequently
been confirmed in more consistent calculations with both valence-space and no-
core methods (cf. Fig. 6).

Multiple studies for $sd$-shell nuclei with a progressively more consistent
perturbative construction of the valence-space interaction followed [134, 135,
136, 204, 141], including first works for the derivation of multi-shell
interactions [138, 144], until the non-perturbative VS-IMSRG and CCEI/SMCC
emerged [139, 140, 142, 143, 179]. The VS-IMSRG, in particular, has been
widely used to compute ground and excited-state energies [142, 143, 129, 225,
226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238], although
theoretical uncertainties stemming from the method still prove challenging
(see Sec. 6).

Figure 9: Deviation from experiment for excited states throughout the $sd$
shell, obtained with _Left:_ the EM 1.8/2.0 interaction without
transformation, _Middle:_ EM 1.8/2.0 transformed with the VS-IMSRG using
ensemble normal ordering, _Right:_ the USDB interaction.

In Fig. 9, we show results from a VS-IMSRG survey of 391 levels in the
$sd$-shell, starting from the EM1.8/2.0 chiral two- plus three-nucleon
interaction [220]. The points indicate the deviation between the computed and
experimental energies for all of these levels, which contribute to the
specified cumulative root-mean-squared deviation between theory and
experiment. In the left panel of Fig. 9, we have simply used the “bare” matrix
elements of the EM1.8/2.0 interaction in the $sd$-shell valence space, while
the center panel shows the results obtained by applying the VS-IMSRG with
ensemble normal ordering, as described in Secs. 2 and 4.1. Since our starting
interaction has been evolved to a low resolution scale, correlations due to
the strong short-range repulsion and the tensor force have largely been
accounted for. Thus, the shell model picture is reasonable: Low-lying nuclear
states are bound and excitation energies are at least of the correct order of
magnitude, with a sizable rms deviation of 1696 keV.

The deviations from experiment are reduced significantly when we use the VS-
IMSRG to decouple the $sd$-shell valence space from other excitations,
accounting for core polarization and other types of long-range, many-body
correlations (see Sec. 2). With an rms deviation of 647 keV, we are not doing
as well as the gold-standard USDB interaction for which the deviation is
merely 220 keV for the selected levels (and only $\sim 130$ keV for all 600+
measured $sd$-shell levels). This is not unexpected: USDB is essentially the
best possible fit to experimental data under the chosen model assumptions,
including the choice of a valence space containing only the
$1s_{1/2},0d_{3/2}$ and $0d_{5/2}$ orbitals, the mass-dependence of the two-
body matrix elements, and the omission of residual three- and higher-body
effective interactions. The accuracy of the VS-IMSRG results is subject to the
uncertainties of the input interaction and the truncation used in the method.
Both can and will be systematically improved in future applications.

###  5.2 The Calcium Region

Soon after the successful application of perturbatively constructed shell
model interactions in the $sd$-shell, first results for the calcium isotopes
followed, including a successful prediction of the two-neutron separation
energies at the sub-shell closure in ${}^{52}\mathrm{Ca}$ [133, 239, 137].
More recently, the masses of ${}^{55-57}\mathrm{Ca}$ were measured at RIKEN,
showing the onset of a flat trend in the separation energies beyond
${}^{54}\mathrm{Ca}$ that would be consistent with the filling of the neutron
$0f_{5/2}$ shell (see Fig. 10). Such a trend had also been found in MR-IMSRG
and Gor’kov Greens Function (GGF) ground-state calculations using chiral
interactions [212], although absolute two-neutron separation values ($S_{2n}$)
could not be determined precisely because of theoretical uncertainties in the
interactions and the many-body methods.

Figure 10: _Left:_ Two-neutron separation energies of neutron-rich calcium
isotopes from recent measurements at RIKEN, compared to VS-IMSRG and MBPT
results obtained with the EM1.8/2.0 interaction, as well as results for
phenomenological interactions. Adapted from [235]. _Right:_ VS-IMSRG (solid
lines) results for three-point energy differences in the calcium isotopes and
neighboring chains, compared to both AME data and new titanium measurements at
TITAN. VS-IMSRG results used the EM1.8/2.0 interaction, while Gor’kov Green’s
Function results for the scandium chain (dashed line) use a different chiral
interaction. Reprinted from [233]. See original references for additional
details.

Recent high-precision mass measurements of the titanium isotopes at TITAN
aimed to shed new light on the evolution of the $N=32$ shell closure [233].
The right panel of Fig. 10 shows three-point energy differences
$\Delta_{2n}\equiv S_{2n}(N,Z)-S_{2n}(N+2,Z)$ extracted from the new data
alongside AME data, in comparison with results from the VS-IMSRG for isotopic
chains in the lower $pf$-shell and GGF calculations for the scandium chain.
While the theoretical $\Delta_{2n}$ compare favorably with experimental data
overall, the strength of the $N=32$ closure is overestimated with increasing
$Z$. This artificial enhancement of shell closures is frequently observed in
calculations with current chiral interactions (see, e.g., [191, 212, 240]),
and might provide important clues toward the refinement of next-generation
forces.

###  5.3 Heavy Nickel and Light Tin

For sufficiently soft interactions, IMSRG and CC calculations for nuclei in
the upper $pf$ and lower $sdg$ shells can be converged [126, 128, 129]. The
limiting factor is a truncation in the three-body matrix elements
$e_{1}+e_{2}+e_{3}\leq E_{3\mathrm{max}}$ where $e=2n+\ell$. Memory
constraints have restricted calculations to $E_{3\mathrm{max}}\leq 18$. The
dimension of the valence space also becomes an issue during the
diagonalization of the effective interaction, but approaches like the Monte
Carlo shell model [241] or importance-truncated configuration interaction (IT-
CI) [242] can be used to tackle this problem.

In Figure 11, we show the evolution of the first excited $2^{+}$ state in
neutron-rich nickel isotopes, which serves as a strong indicator for
(sub-)shell closures. The jump in the $2^{+}$ energy at ${}^{78}\mathrm{Ni}$
suggests that this nucleus is indeed doubly magic. The VS-IMSRG reproduces the
available experimental data [243] well, and the energies are insensitive under
(admittedly small) variations of the interaction’s resolution scales or low-
energy constants — see [220] for more details on these Hamiltonians.

Recently, Hagen _et al._ also computed the $2^{+}$ [128] energies of
${}^{78,80}\mathrm{Ni}$ using the Equation-of-Motion Coupled Cluster approach.
In Fig. 11, we include their excitation energies from the so-called EOM-
CCSD(T) method, which are about 1 MeV lower than the VS-IMSRG results with the
corresponding interactions. This difference can be traced back to the effects
of triples (i.e., 3p3h) correlations and continuum effects that are currently
not included in the VS-IMSRG.

Figure 11: Energies of the first excited $2^{+}$ states in Ni isotopes from
VS-IMSRG and Coupled Cluster calculations including triples corrections [128],
using EM1.8/2.0 and other chiral two- plus three-nucleon interactions as input
(see [220] for details).

Moving to even heavier nuclei, the structure of the lightest tin isotopes was
the subject of a recent joint EOM-CC and VS-IMSRG study [129]. Figure 12 shows
results for the energy gap between the two lowest-lying states in light odd-
mass tin isotopes and ${}^{105}\mathrm{Te}$. The no-core EOM-CC and the VS-
IMSRG results for ${}^{101}\mathrm{Sn}$ are consistent, and the VS-IMSRG
produces a near-degeneracy of the $J^{\pi}=5/2^{+}$ and $J^{\pi}=7/2^{+}$
states that is compatible with experiment, the systematic uncertainties of the
method must be properly quantified — and, most likely, reduced — before one
can make spin assignments with confidence.

Figure 12: VS-IMSRG and EOM-CC results for the ground and first excited states
of odd-mass tin isotopes and ${}^{105}\mathrm{Te}$, using chiral NN+3N
interactions [220]. The right panel shows the convergence of the states as a
function of the model space truncation in the IT-CI diagonalization. Figure
adapted from Ref. [129].

##  6 Current challenges

While great strides have been made in deriving effective interactions for the
shell model, challenges remain. Here we focus on two in particular, and
analyze them from the perspective of the VS-IMSRG.

First, electric quadrupole (E2) observables which are sensitive to low-lying
collective excitations and which historically have been treated
phenomenologically by introducing effective charges, are not captured well
with present techniques. This can qualitatively be understood in the context
of the cluster expansion discussed at the beginning of Sec. 2. It is precisely
the low-lying _collective_ modes that are expected to violate the cluster
hierarchy upon which the IMSRG relies.

Second, several regions of the nuclear chart—e.g. the “islands of inversion”
[244, 245], or the charge radii of the calcium isotopes [246, 247]—display
features which suggest that a naive valence space of a single major harmonic
oscillator shell is not an appropriate first-order description. However, the
derivation of effective interactions for non-standard valence spaces leads to
difficulties related to the well-known intruder-state problem.

###  6.1 Electromagnetic transitions

The first attempt at a microscopic treatment of electric quadrupole ($E2$)
observables was the work of Horie and Arima [248], investigating the role of
configuration mixing on quadrupole moments. A series of investigations by
Seigel and Zamick [249, 250, 251] demonstrated the importance of terms beyond
first order in perturbation theory. Specifically, they investigated the impact
of Tamm-Dancoff (TDA) and random-phase approximation (RPA) graphs to the
effective charge, with the physical interpretation that the effective charge
comes largely from a coupling to the giant quadrupole resonance. A subsequent
calculation by Kirson [44] indicated that a self-consistent treatment
including screening effects essentially canceled the effect obtained with RPA.
For a discussion, see Ref. [27].

An important development came with the application of the OLS approach to an
effective interaction for 6Li in the $p$-shell [162], where the resulting
effective $E2$ operator could be reasonably approximated by the use of
effective charges. This approach was investigated in more detail over a decade
later [164], showing that a nonperturbative treatment could produce the
collective effects of $E2$ observables. Of course, it is not always clear what
lessons learned in the $p$ shell carry over to heavier masses.

As discussed in Sec. 2.3.2, the Magnus formulation of the IMSRG provides a
straightforward way to construct effective valence space operators for general
observables. All operators are consistently transformed according to

$\mathcal{O}_{\mathrm{eff}}=e^{\Omega}\mathcal{O}e^{-\Omega}=\mathcal{O}+[\Omega,\mathcal{O}]+\tfrac{1}{2}\left[\Omega,[\Omega,\mathcal{O}]\right]+\ldots$
(77)

A first application of this approach was to electromagnetic transitions in
light and medium mass nuclei [123], where it was found that the observables
were well-converged with respect to the model space truncation (i.e. frequency
and number of major shells included in the initial harmonic oscillator basis).
However, the computed values for collective observables like magnetic moments
or electric quadrupole and octupole transitions were substantially smaller
than experimental data.

The possible explanations for this discrepancy are that either that the
truncation of Eq. (77) to two-body operators is insufficient to capture this
type of collectivity, or else the input chiral interactions are deficient in
some way. Most likely, both are in effect to some degree. The interaction used
in Ref. [123] is known to underpredict charge radii in these same nuclei
[240]. Given that the electric quadrupole operator is proportional to $r^{2}$,
where $r$ is the point proton radius, and that the transition strengths
$B(E2)$ go as $r^{4}$, one would naturally expect some underestimation of the
quadrupole strength. However, as demonstrated in Fig. 13, this cannot be the
whole story.

Figure 13: Electric quadrupole transition matrix element $\langle
0^{+}\|E2\|2^{+}\rangle$ in 14C, 22O, and 32S computed using the VS-IMSRG with
several choices of input chiral interaction. Also shown are the energy of the
$2^{+}_{1}$ state and the point proton radius squared. Experimental radii are
from [252], energies and transition matrix elements are from [243].

The point proton radius squared, indicated $R^{2}_{pp}$ in Fig. 13, is
underpredicted at approximately the same level in 14C and 32S. On the other
hand, while the $E2$ strength for 14C is reasonably reproduced, in 32S it is
underpredicted by $\sim$25-50%, and the strength in 22O is underpredicted by
$\sim$65%. Clearly, the underprediction of $E2$ strength in 32S cannot be
explained solely by the radius deficiencies.

Supporing this interpretation, unpublished calculations in a small space where
exact diagonalizations are possible show unambiguously that the IMSRG(2)
truncation misses a significant fraction of the $E2$ strength, and that
capturing the full strength requires inclusion of correlated many-particle
many-hole excitations. In addition, symmetry-adapted no-core shell model
calculations of $sd$ shell nuclei yielded significantly larger $E2$ strength
using the same starting interaction [231].

It is interesting to compare the IMSRG results with the previously-mentioned
approach of Siegel and Zamick [251]. They considered three levels of
approximation for the $E2$ operator, first-order core-polarization, TDA, and
RPA. Typical diagrams contributing in these approximations are shown in Fig.
14. The effective operator generated by the VS-IMSRG via Eq. (77) contains TDA
and RPA graphs to all orders. The VS-IMSRG also sums higher-order diagrams
like Fig. 14(d), but because of the truncation to two-body operators, certain
types of diagrams are undercounted or missing altogether (see Refs. [145, 174,
175] for more details).

(a) Core pol.

(b) TDA

(c) RPA

(d) IMSRG

Figure 14: Exaples of diagrams contributing to the one-body part of the effective $E2$ operator. The X represents the bare operator. Table 1: Neutron effective charges for the $sd$ shell obtained with first-order core polarization, including TDA and RPA graphs to all orders, and from IMSRG. These results are obtained with the EM1.8/2.0 interaction [220], in a Hartree-Fock basis constructed from an oscillator basis with $e_{max}=10$, $\hbar\omega=16\,\text{MeV}$. | $d_{5/2}d_{5/2}$ | $d_{3/2}d_{3/2}$ | $d_{5/2}d_{3/2}$ | $d_{5/2}s_{1/2}$ | $d_{3/2}s_{1/2}$
---|---|---|---|---|---
Core. pol. | 0.110 | 0.035 | 0.064 | 0.034 | 0.026
TDA | 0.121 | 0.037 | 0.062 | 0.040 | 0.031
RPA | 0.119 | 0.037 | 0.061 | 0.038 | 0.030
IMSRG | 0.202 | 0.098 | 0.222 | 0.163 | 0.093

Table 1 presents the effective charge for a neutron in the $sd$ shell in these
various levels of approximation. The orbit-dependent effective charge is
obtained as [251]

$e_{ab}=\frac{\langle a\|\mathcal{O}^{E2}_{\mathrm{eff}}\|b\rangle}{\langle\pi
a\|\mathcal{O}^{E2}_{\mathrm{bare}}\|\pi b\rangle}$ (78)

where in the denominator, we take the matrix element for the corresponding
proton orbit. We work in a Hartree-Fock basis constructed from an oscillator
basis with frequency $\hbar\omega=16\,\text{MeV}$ and $e_{max}=10$. In this
basis, we obtain a bare proton matrix element $\langle\pi
d_{5}\|\mathcal{O}^{E2}_{\mathrm{bare}}\|\pi
d_{5}\rangle=-9.03e\,\text{fm}^{2}$, and we see that in order to reproduce the
experimental quadrupole moment of 17O ($Q=-2.56(2)e\,\text{fm}$ [253]), we
require an effective neutron charge of161616With our definition of the $E2$
operator, $Q=\sqrt{16\pi/5}\langle J,M=J|\mathcal{O}^{E2}|J,M=J\rangle$
$e_{n}\approx 0.37$. Likewise, the bare proton matrix element $\langle\pi
d_{5}\|\mathcal{O}^{E2}_{\mathrm{eff}}\|\pi
s_{1}\rangle=-9.29e\,\text{fm}^{2}$, and so to reproduce the experimental
transition strength
$B(E2;\tfrac{1}{2}^{+}\rightarrow\tfrac{5}{2}^{+})=6.21(8)e^{2}\,\text{fm}^{4}$
we require an effective neutron charge $e_{n}\approx 0.38$. We find that while
the IMSRG generates a larger neutron effective charge than the other methods,
the result is still well below the experimental value.

These IMSRG effective charges are essentially the same as those found in a
previous study [123] with a different chiral interaction. That study also
found proton effective charges close to 1, i.e. with almost no
renormalization. As discussed in Ref. [27], this can potentially be understood
by considering that in order to “dress” a valence nucleon, that nucleon must
excite a proton out of the core. A valence neutron can do this through the
$T=0$ channel, while a valence proton must act in the weaker $T=1$ channel.

###  6.2 The Intruder-state Problem

Over the last few decades, experimental investigations of nuclei far from
stability have revealed the existence of several “islands of inversion”, where
nuclei near traditional shell closures have ground states that indicate
significant deformation or correlated particle-hole excitations out of the
closed shell [254, 244, 245]. The classic examples are 31Na and 32Mg, both
with $N=20$. 31Na has a ground state spin-parity of $\tfrac{3}{2}^{+}$ , while
shell model calculations predicted $\tfrac{5}{2}^{+}$, and 32Mg has a $2^{+}$
excitation energy of 885 keV, far lower than expected for a closed neutron
shell. Both have greater binding energies than predicted in the shell model.
If these correlations are sufficiently important, then it is possible that the
ground state will not be among the subset of eigenstates reproduced in the
valence space diagonalization. Indeed, Watt et al. [255] found that by
explicitly allowing neutron excitations out of the standard $sd$ shell and
into the $f_{7/2}$ shell, the discrepancies for 31Na and 32Mg could be
understood. Of course, even if the correlated ground state is formally among
the valence states, it is likely that incorporating the correlated excitations
into an effective Hamiltonian would require large many-body forces.

It is therefore desirable to be able to produce an effective interaction for a
valence space that spans more than one major shell, such as the $sd-fp$ space,
and indeed phenomenological interactions for such a space have been successful
at describing the island of inversion effects [245], as well as the charge
radii of the calcium isotopes [246]. Unfortunately, in deriving such an
interaction from first principles, one runs into the well-known intruder-state
problem which we discuss below. In fact, attempts to include effects of the
continuum—essential for studies near the driplines—suffer from the same
problem [256]. Understanding and solving this problem, particularly in the
context of a nonperturbative approach, is clearly of great interest.

####  6.2.1 The intruder-state problem in perturbation theory

As demonstrated by Schucan and Weidenmüller [46, 47], there are serious
reasons to doubt the convergence of the perturbative expansion for the
effective interaction. To illustrate this, we split up the Hamiltonian as
before into a zero-order piece and a perturbation, with a dimensionless power-
counting parameter $x$

$H(x)=H_{0}+xV$ (79)

such that $H(0)$ is the zero-order Hamiltonian and $H(1)$ is the full
Hamiltonian. The perturbative expansion of $H_{\mathrm{eff}}$ can therefore be
seen as a Taylor expansion about $x=0$ evaluated at $x=1$. The trouble arises
if one of the states belonging to the $Q$ space has an energy lower than one
of the $P$ space states. Such a state is called an “intruder state”. Assuming
the $P$ states are all at lower energy than the $Q$ states at $x=0$, this
implies a level crossing for some $x\in[0,1]$. The value of $x$ at which such
a crossing occurs (even if it is an avoided crossing) corresponds to a branch
point which places an upper limit on the radius of convergence of the
effective Hamiltonian [47].

Unfortunately such level crossings are the rule, not the exception. Moreover,
if the zeroth-order levels in the valence space are non-degenerate—e.g. if a
Hartree-Fock basis is used—then as more particles are added to the valence
space, the energy of the highest $P$ space can quickly become higher than the
energy of the lowest $Q$ space, even without the residual interaction.

####  6.2.2 The intruder-state problem in IMSRG

There are a few reasons why one might expect the IMSRG to avoid the intruder-
state problem. First, it is formally a non-perturbative method, so the above
argument does not directly apply. Second, because it is formulated in Fock
space, the energies of the $A$-body system do not enter into any energy
denominators, and so one would not naively expect divergences due to crossings
in the $A$-body system.

Unfortunately, the IMSRG suffers from a related, but distinct intruder-state
problem. An illustrative example of the type of behavior encountered is shown
in Fig. 15. Here we aim to decouple a valence space consisting of the $p$ and
$sd$ major shells from a large space constructed from 7 major harmonic
oscillator shells ($e_{max}=6$). We use an 16O Hartree-Fock reference state,
which is indicated schematically in Fig. 15. Also shown in Fig. 15 are the
zero-body term $E_{0}(s)$, the norm of the generator $\|\eta(s)\|$ and the
norm of the Magnus operator $\|\Omega(s)\|$ as a function of the flow
parameter $s$. As usual, we do this in two steps, first decoupling excitations
out of the core (4He in this case), followed by a decoupling of the valence
space (cf. Sec. 2.3). The core decoupling is achieved at $s\approx 12$. At
this point, we observe a jump in $\|\eta(s)\|$ because our new definition of
‘off-diagonal” now includes many more matrix elements. In a well-behaved
calculation, these terms would then be suppressed by the IMSRG flow. Indeed,
the size of $\eta$ initially decreases, but it soon begin to grow again, and
the caclulation fails to converge. We also observe that the flow of the zero-
body term $E_{0}$ turns around and diverges, and the Magnus operator $\Omega$
grows indefinitely. At some point, $\Omega$ grows beyond the radius of
convergence of the BCH expansion. As a result, no effective interaction is
obtained.

We mention in passing that there has been some success using the IMSRG to
decouple valence spaces other than those defined by a single major harmonic
oscillator shell, so long as they are reasonably well separated by a shell
gap. These spaces correspond to the “extruded-intruded” spaces described by
the Strasbourg group [18], where due to the spin-orbit potential the orbit
with the largest-$j$ orbit drops out (is “extruded”) and the largest-$j$ orbit
from the next shell up comes down (it “intrudes”). An example is the space
consisting of the orbits $1p_{3/2}$, $1p_{1/2}$, $0f_{5/2}$, $0g_{9/2}$. This
space (for neutrons), was used to treat heavy chromium isotopes [236].
However, the results obtained there suggested that this space was not
sufficient to describe the ground states of those isotopes.

Figure 15: _Left:_ Decoupling of the $psd$ valence space using an 16O
reference, shown with the Hartree-Fock single-particle spectrum. _Right:_ The
zero-body part of the flowing Hamiltonian, the norm of the generator $\eta$
and the norm of the Magnus operator $\Omega$ as a function of the flow
parameter $s$. At $s\sim 12$, the core is decoupled and the decoupling of the
valence space begins.

The connection between intruders and failed convergence of the IMSRG flow can
be understood schematically by considering the flow equation formulation
$\tfrac{d}{ds}H=[\eta,H]$. Imagine we have two levels labeled $p$ and $q$,
both with degeneracy greater than 2, and with single particle energies
$\epsilon_{q}>\epsilon_{p}$. We intend to decouple the $q$ orbit from the
Hilbert space, which means suppressing terms like
$V_{qqpp}a^{\dagger}_{q}a^{\dagger}_{q}a_{p}a_{p}$ which excite particles from
the $p$ level to the $q$ level. The flow equation for $V_{qqpp}$ is,
schematically,

$\frac{d}{ds}V_{qqpp}\sim
V_{qqpp}\frac{2\epsilon_{p}-2\epsilon_{q}+V_{pppp}-V_{qqqq}}{2\epsilon_{q}-2\epsilon_{p}}+\ldots$
(80)

If the one-body terms dominate the right hand side of (80), then
$\tfrac{d}{ds}V_{qqpp}\sim-V_{qqpp}$ and the off-diagonal term is suppressed
exponentially. On the other hand, if the interaction terms $V$ are larger than
the one-body terms, and of opposite sign, then $V_{qqpp}$ will be
exponentially enhanced. This can be achieved if $V_{qqqq}$ is negative
(attractive) and $V_{pppp}$ is positive (repulsive), and a positive numerator
in (80) corresponds to an inversion of the states $|pp\rangle$ and
$|qq\rangle$.

In fact, terms like $V_{pppp}$ and $V_{qqqq}$ can be included in the
denominator by repartitioning the Hamiltonian so that the diagonal (i.e.
bra=ket) parts of $V$ are included in $H_{0}$, avoiding the exponential
growth. However, intruders can also be driven by collective effects which
cannot be tamed by a straightforward repartitioning. Consider the case where
we have multiple included and excluded levels $p,p^{\prime}\ldots
q,q^{\prime}\ldots$ In this case, we should also consider contributions like

$\frac{d}{ds}V_{qqpp}\sim\frac{V_{qqp^{\prime}p^{\prime}}}{2\epsilon_{q}-2\epsilon_{p^{\prime}}}V_{p^{\prime}p^{\prime}pp}-V_{qqq^{\prime}q^{\prime}}\frac{V_{q^{\prime}q^{\prime}pp}}{2\epsilon_{q^{\prime}}-2\epsilon_{p}}+\ldots$
(81)

If there are many such terms involving
$V_{p^{\prime\prime}p^{\prime\prime}pp}$, etc. and these terms add coherently,
with the $V_{qqq^{\prime\prime}q^{\prime\prime}}$ type terms having opposite
sign, they can compete with the contributions in (80) and potentially lead to
growth of the off-diagonal terms. Such a situation will also lead to a
crossing of collective levels. Clearly, this situation and the previous one
will be exacerbated by the small energy denominators which occur in multi-
shell valence spaces.

There is another way in which intruders can cause trouble, and this is by
spoiling the cluster hierarchy. We illustrate this with a toy system in the
next section.

####  6.2.3 Toy model for the intruder-state problem

To illustrate how intruders and level crossings can lead to large induced
many-body terms, we consider the problem of three kinds of fermion—which could
be, say, spin-up neutron, spin-down neutron, spin-up proton— living in a
three-level Hilbert space. We require three particles because we wish to
monitor induced three-body forces.

The initial Hamiltonian is $H(x)=H_{0}+xV$ where

$H_{0}\equiv\sum_{i}\epsilon_{i}a^{\dagger}_{i}a_{i}\hskip 10.00002pt,\hskip