††This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-
AC05-00OR22725 with the U.S. Department of Energy. The United States
Government retains and the publisher, by accepting the article for
publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for United
States Government purposes. The Department of Energy will provide public
access to these results of federally sponsored research in accordance with the
DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

# IMSRG with flowing 3 body operators, and approximations thereof

S. R. Stroberg<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Notre Dame, Notre Dame, IN, 46556, USA  Physics Division,
Argonne National Laboratory, Lemont, IL, 60439, USA  Department of Physics,
University of Washington, Seattle, WA 98195, USA  T. D. Morris  Quantum
Information Science Section, Oak Ridge National Laboratory, Oak Ridge, TN
37831, USA  B. C. He  Department of Physics and Astronomy, University of Notre
Dame, Notre Dame, IN, 46556, USA

###### Abstract

We explore the impact of retaining three-body operators within the in-medium
similarity renormalization group (IMSRG), as well as various approximations
schemes. After studying two toy problems, idential fermions with a contact
interaction and the Lipkin-Meshkov-Glick model, we employ the valence-space
formulation of the IMSRG to investigate the even-$A$ carbon isotopes with a
chiral two-body potential. We find that retaining only those commutators
expressions that scale as $N^{7}$ provides an excellent approximation of the
full three-body treatment.

††preprint: APS/123-QED

##  I Introduction

There has been considerable progress in the past few years in ab initio many-
body methods for atomic nuclei [1], with an explosion in the number of
available methods, as well as in the systems amenable to ab initio treatment.
At present, nuclei from $A=2$ up to $A\sim 208$ can be accessed [2, 3, 4],
including many open-shell nuclei [5, 6, 7, 8, 9]. There has simultaneously
been great progress in quantifying uncertainties due to the truncation of
effective field theory (EFT) to a finite order, especially the use of
emulators for the solution of the many-body problem [10, 11, 12, 13, 14]. One
of the most important remaining open issues in ab initio nuclear theory is a
robust and reliable assessment of the error due to approximations made in
solving the many-body problem.

In this paper, we focus on the in-medium similarity renormalization group
(IMSRG) approach [15, 16, 17, 18, 19], and work towards understanding how the
adopted truncation scheme manifests as errors in the predicted observables.
The paper is organized as follows. In section II, we summarize the IMSRG
method and discuss the truncation scheme. In section III we discuss approaches
to approximate a full treatment of three-body operators within the IMSRG. In
section IV we explore two toy models to develop intuition of the IMSRG
truncation scheme, before applying the machinery to a more realistic
calculation of the carbon isotopes. Finally, in section V we briefly discuss
the computational cost of our implementation of the IMSRG.

##  II In-medium similarity renormalization group

In this section, we briefly recapitulate the formalism of the in-medium
similarity renormalization group, which has previously been summarized in
several review papers, including [17, 19].

###  II.1 Basic ingredients

The basic approach of the SRG is to perform a unitary transformation on the
Hamiltonian

$H(s)\equiv U(s)HU^{\dagger}(s),$ (1)

in order to bring it to a form which is more amenable to solving the many-body
Schrödinger equation. The dependence of the unitary transformation $U(s)$ on
the flow parameter $s$ is specified by the generator $\eta$

$\frac{d}{ds}U(s)=\eta(s)U(s).$ (2)

To ensure that $U(s)$ is unitary, we require that $\eta$ is anti-hermitian,
$\eta^{\dagger}=-\eta$. This leads to the SRG flow equation

$\frac{d}{ds}H(s)=[\eta(s),H(s)].$ (3)

Integration of the ordinary differential equation (3) yields $H(s)$. Given
that the only requirement on $\eta$ is that it is anti-hermitian, there is
some art in choosing a generator which efficiently leads to an improved form
of the Hamiltonian. Here, we use the White generator [20], or its arctangent
variant which we write schematically as

$\eta^{\rm Wh}=\frac{H^{\rm od}}{\Delta}\hskip 10.00002pt,\hskip
20.00003pt\eta^{\rm atan}=\frac{1}{2}{\rm atan}\left(\frac{2H^{\rm
od}}{\Delta}\right).$ (4)

In (4) $H^{\rm od}$ is the suitably-defined off-diagonal part of the
Hamiltonian (the part which should be suppressed by the RG evolution).
$\Delta$ is an energy denominator, typically given by Epstein-Nesbet or
Møller-Plesset partitioning. We employ a superoperator notation; the division
and arctangent should be understood as acting elementwise. For a more explicit
description, see e.g. [21].

###  II.2 Magnus formulation

Rather than directly integrating (3), it is often convenient both numerically
and formally to use the Magnus formulation of the SRG [22]. In this case, we
express the unitary transformation as an exponentiated anti-hermitian operator
$\Omega$

$U(s)=e^{\Omega(s)}.$ (5)

We may derive a flow equation for $\Omega(s)$ by using the Baker-Campbell-
Hausdorff formula

$\frac{d}{ds}\Omega(s)=\sum_{k=0}^{\infty}\frac{B_{k}}{k!}[\Omega(s),\eta(s)]^{(k)}$
(6)

where $B_{k}$ is the $k$th Bernoulli number and $[\eta,\Omega]^{(k)}$
indicates a $k$-fold nested commutator

$[\Omega,\eta]^{(k)}=\begin{cases}[\Omega,[\Omega,\eta]^{(k-1)}]&k>0\\\
\eta&k=0\end{cases}.$ (7)

The flowing Hamiltonian is then given by

$H(s)=e^{\Omega(s)}H(0)e^{-\Omega(s)}=\sum_{k=0}^{\infty}\frac{1}{k!}[\Omega(s),H(0)]^{(k)}.$
(8)

In order to perform a calculation, we must choose a representation for our
operators. Typically we choose a Fock space representation, writing operators
in terms of strings of creation and annihilation operators. For example, a
particle-number conserving operator $O$ would be written as

$\displaystyle O=O_{0}$
$\displaystyle+\sum_{ij}O_{ij}\\{a^{\dagger}_{i}a_{j}\\}+\frac{1}{(2!)^{2}}\sum_{ijkl}O_{ijkl}\\{a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}\\}$
(9)
$\displaystyle+\frac{1}{(3!)^{2}}\sum_{ijklmn}O_{ijklmn}\\{a^{\dagger}_{i}a^{\dagger}_{j}a^{\dagger}_{i}a_{m}a_{n}a_{l}\\}+\ldots$

where the ellipses indicate terms involving more than three creation and three
annihilation operators. The braces $\\{\cdot\\}$ indicate that the the string
of creation and annihilation operators are normal ordered with respect to some
reference state $|\Phi\rangle$ (which could possibly be the true vacuum
$|0\rangle$) so that $\langle\Phi|\\{a^{\dagger}\ldots a\\}|\Phi\rangle=0$.
The coefficients $O_{ijkl},O_{ijklmn}$ are antisymmetrized matrix elements so
that, e.g., $O_{ijkl}=\langle ij|O|kl\rangle$, where
$|kl\rangle=a^{\dagger}_{k}a^{\dagger}_{l}|0\rangle$.

Occasionally, we find that the norm $\|\Omega\|$ becomes large during the
IMSRG flow and evaluating (6) and (8) requires the evaluation of many nested
commutators before convergence is reached. In this case it becomes
computationally advantageous to split up the transformation so that (5)
becomes

$U(s)=U(s-s_{1})U(s_{1})=e^{\Omega(s-s_{1})}e^{\Omega(s_{1})}.$ (10)

We then save the intermediate Hamiltonian $H(s_{1})$ and evaluating
$H(s)=e^{\Omega(s-s1)}H(s_{1})e^{-\Omega(s-s1)}$ converges more rapidly. If
$\Omega(s-s_{1})$ becomes large, we can repeat the procedure so that after $n$
splits,

$U(s)=e^{\Omega(s-s_{n})}\ldots e^{\Omega(s_{2}-s_{1})}e^{\Omega(s_{1})}.$
(11)

If we are only interested in energies, we only need to retain the most recent
$H$ and $\Omega$ operators. In the limit that we split after each
infinitessimal step $ds$, the procedure becomes equivalent to directly
integrating the flow equation (3).

For certain approximation schemes, such as the IMSRG(2*) approximation
described below, this splitting leads to an enhanced error due to missing
cross terms. It is then useful to adopt a strategy which we call “hunter-
gatherer”, in which we split the transformation into two parts

$U(s)=e^{\Omega_{G}(s)}e^{\Omega_{H}(s)}$ (12)

where $\Omega_{H}$ is the “hunter” and $\Omega_{G}$ is the “gatherer”.
Whenever the norm $\|\Omega_{H}\|$ gets beyond a (typically small) threshold,
we update the gatherer via the Baker-Campbell-Hausdorff formula

$\Omega_{G}\to\sum_{k=0}^{\infty}\frac{B_{k}}{k!}[\Omega_{G},\Omega_{H}]^{(k)}.$
(13)

This approach retains the benefits of the aforementioned splitting strategy,
without losing the cross terms in the IMSRG(2*) approximaiton.

###  II.3 Truncation of the SRG flow equations and the NO$n$B approximation

From (9) combined with either (3) or (6) and (8), we can see that the main
computational task required for an IMSRG calculation is the evaluation of
commutators of Fock-space operators. Even if $\eta$ and $H$ are initially two-
body operators, integrating (3), or evaluating the nested commutators in (7)
and (8) will induce three-body and higher-body operators. To make the
calculation feasible, all operators, including intermediate commutator
expressions, are truncated at the two-body level, yielding the IMSRG(2)
approximation111In free-space SRG calculations, due to the translation-
invariance of the reference $|\Phi\rangle=|0\rangle$, three-body operators can
be managed more efficiently by eliminating the center of mass coordinate and
working in Jacobi coordinates [23, 24])..

The IMSRG(2) approximation can be viewed as repeatedly making the normal-
ordered two-body (NO2B) approximation, which is also used with other methods,
and which is very accurate beyond $A\sim 16$ [25, 26, 27]. For a single-
determinant reference $|\Phi\rangle$, the normal-ordered $n$-body (NO$n$B)
approximation becomes exact when $n$ is equal to the maximum number of quasi-
particles (particles plus holes) present in the exact wave function. For
example, a one-particle-one-hole (two quasi-particle) configuration receives
contributions from a normal-ordered two-body operator; but a normal ordered
three-body operator will not contribute:
$\langle\Phi^{a}_{i}|\\{a^{\dagger}a^{\dagger}a^{\dagger}aaa\\}|\Phi^{b}_{j}\rangle=0$.
Likewise, a two-particle-two-hole configuration receives contributions from a
four-body operator, but not a five-body operator.

For an $A$-body system with a single-determinant reference with $A$ particles,
at most $A$-particle-$A$-hole configurations can be present in the wave
function, and so the NO$n$B approximation is exact for $n\geq 2A$ [28]. With a
vacuum reference, at most $A$ quasiparticles are present, and so the NO$n$B
approximation is exact for $n\geq A$. We are interested in treating systems
with $A\gg 3$, and going beyond the NO3B approximation is generally
prohibitively expensive, so the condition $n\ll A$ is inevitable.

The NO$n$B approximation will also be exact if there are no $m$-body operators
with $m>n$; for example, working at low order in chiral effective field theory
[29, 30] there are no operators with $m>3$, so the NO3B approximation is
exact. However, SRG evolution will induce operators with $m>3$ so IMSRG(3) is
no longer exact. The error due to neglecting $m$-body operators will depend on
the relative importance of $m$-quasiparticle configurations in the wave
function. For the SRG with a vacuum reference, all configurations have $A$
quasiparticles, while for IMSRG (with a judicious choice of reference) the
wave function should be dominated with $m$-quasiparticle configurations with
$m\ll A$.

Focusing for concreteness on the NO2B approximation of a three-body operator,
we schematically write the expectation value of a string of
creation/annihilation operators in the exact wave function $|\Psi\rangle$,
with reference $|\Phi\rangle$, as

$\displaystyle\langle\Psi|a^{\dagger}a^{\dagger}a^{\dagger}aaa|\Psi\rangle\sim{}$
$\displaystyle\langle\Phi|a^{\dagger}a^{\dagger}a^{\dagger}aaa|\Phi\rangle$
(14) $\displaystyle+$
$\displaystyle\langle\Phi|a^{\dagger}a^{\dagger}aa|\Phi\rangle~{}\langle\Psi|\\{a^{\dagger}a\\}|\Psi\rangle$
$\displaystyle+$
$\displaystyle\langle\Phi|a^{\dagger}a|\Phi\rangle~{}\langle\Psi|\\{a^{\dagger}a^{\dagger}aa\\}|\Psi\rangle$
$\displaystyle+$
$\displaystyle\langle\Psi|\\{a^{\dagger}a^{\dagger}a^{\dagger}aaa\\}|\Psi\rangle.$

The first term in (14) will scale as $A^{3}$. The second term will scale as
$A^{2}N_{q}$, where $N_{q}$ is the number of quasiparticles 222Here, we assume
that the size of $\langle\Psi|\\{a^{\dagger}a\\}|\Psi\rangle$ is reasonably
approximated by the diagonal terms $N_{q}=\sum_{p}\langle
p^{\dagger}p\rangle+\sum_{h}\langle hh^{\dagger}\rangle$. In nuclear structure
applications, many of the off-diagonal terms will be suppressed by selection
rules on the quantum numbers $j,\pi,t_{z}$.. The third term will go as
$AN_{q}^{2}$, and the fourth as $N_{q}^{3}$, so we see that we have an
organization in powers of $N_{q}/A$ [31, 32]. The approximation is improved by
the IMSRG flow, which suppresses excitations so that $N_{q}(s)\to 0$ as
$s\to\infty$. This argument is complicated somewhat by the fact that the
three-body operators induced by the IMSRG flow are not uniformly distributed
over all orbits, but instead are concentrated near the Fermi surface.
Nevertheless, so long as the typical number of quasiparticles is smaller than
the number of particles feeling the induced forces, we expect the effect of
many-body operators to be suppressed.

Figure 1: Hugenholtz skeleton diagrams indicating the commutator topologies
present in IMSRG(3) which are omitted in the IMSRG(2) approximation. The
circle and square represent the two operators entering into the commutator.

##  III Approximations to full IMSRG(3)

The full commutator expressions necessary for IMSRG(3) in the uncoupled and
$J$-coupled representations are presented in Appendices A and B respectively
(see also refs [33, 34, 17, 35]). The corresponding Hugenholtz skeleton
diagrams which enter at the IMSRG(3) level are presented in Fig. 1. The naive
computational scaling to evaluate a particular term can be read off of the
diagram by counting the number of fermion lines. For example, the first
diagram in Fig. 1, labeled333Here we use the shorthand notation $[a,b]_{c}$
for the $c$-body piece of the commutator of an $a$-body operator with a
$b$-body operator. $[3,3]_{0}$, has six fermion lines and thus scales as
$N^{6}$, where $N$ is the number of single-particle states in the basis. The
final diagram in Fig. 1 involves nine fermion lines and thus scales as
$N^{9}$. This is too expensive for realistic applications (see section V), and
so we shall explore some approximation schemes which can render the
calculation more tractable.

###  III.1 Goose-tank diagrams and the IMSRG(2∗) approximation

In refs. [34, 17], a perturbative analysis of the IMSRG demonstrated that the
IMSRG(2) ground state energy is exact through third order in MBPT, while some
fourth-order and higher-order terms are missed. Specifically, all fourth-order
diagrams with three-particle-three-hole intermediate states (triples) are
missed, while a subset of diagrams with four-particle-four-hole intermediate
states (asymmetric quadruples, shown in Fig. 2) are undercounted by a factor
$1/2$. Coupled cluster with singles and doubles (CCSD) also misses the fourth-
order triples, but correctly includes the quadruples [34].

Figure 2: Hugenholtz skeleton diagrams indicating the fourth-order quadruples
which are undercounted in the IMSRG(2) (cf Fig 5.5 of [36]).

Within the IMSRG, the energy diagrams in Fig. 2 are obtained444The commutator
of the third-order generator with the first order vertex yields an identical
expression[34, 17]. from the commutator of the third order two-body vertex
$\Gamma^{[3]}$ with the first-order generator $\eta^{[1]}$. The missing half
of the quadruples diagrams are due to the contributions to $\Gamma^{[3]}$
which involve an intermediate three-body operator. These are proportional to
$[\eta,[\eta,\Gamma]_{\rm 3b}]_{\rm 2b}$, as indicated in Figs. 3(a) and 3(b),
where $\eta$ and $\Gamma$ are evaluated at $s=0$. We refer to these diagrams
as “goose-tank” diagram because of their appearance when rotated 90 degrees.
We note that their importance was also identified by Evangelista and
collaborators in the context of quantum chemistry [37, 38, 39].

The half of the asymmetric quadruples which are included in the IMSRG(2) arise
from contributions to $\Gamma^{[3]}$ proportional to $[\eta,[\eta,\Gamma]_{\rm
1b}]_{\rm 2b}$, involving an intermediate one-body operator, as indicated in
Figs. 3(c) and 3(d). Evaluation of the diagrams in Fig. 3 reveals that 3(c) is
equal to 3(a), while 3(d) is only approximately equal to 3(b). Nevertheless,
3(d) and 3(b) give identical contributions to the fourth-order energy when
further contracted with the generator $\eta$ to yield a zero-body operator.
The full fourth-order quadruples can therefore be restored by modifying the
IMSRG(2) to the so-called IMSRG(2∗) scheme. In the standard formulation, we
modify the flow equation to be

$\frac{d}{ds}H(s)=[\eta(s),H(s)+\chi(s)]$ (15)

where $\chi(s)$ is an auxiliary one-body operator which obeys the flow
equation

$\frac{d}{ds}\chi_{ij}=(n_{i}n_{j}+\bar{n}_{i}\bar{n}_{j})[\eta_{\rm
2b},H_{\rm 2b}]_{ij}.$ (16)

In the Magnus formulation, the nested commutators used in (8) are modified to

$[\Omega,H]^{(k+1)}=\begin{cases}[\Omega,[\Omega,H]^{(k)}]&k=0\\\
[\Omega,[\Omega,H]^{(k)}+\chi^{(k)}]&k>0,\end{cases}$ (17)

and the one-body operator $\chi^{(k)}$ is defined as

$\chi^{(k)}_{ij}=(n_{i}n_{j}+\bar{n}_{i}\bar{n}_{j})[\Omega_{\rm
2b},[\Omega,H]^{(k-1)}_{\rm 2b}]_{ij}.$ (18)

As noted in section II.2, the procedure is made more complicated if we split
the transformation into multiple steps, and so it is best to use the hunter-
gatherer scheme.

Figure 3: (a) and (b): Diagrammatic illustration of a contributions to
$[\eta,[\eta,H]_{3}]_{2}$, which we denote “goose-tank” diagrams. (c) Diagram
arising from $[\eta,[\eta,H]_{1}]_{2}$ which gives a contribution identical to
that of (a). (d) Diagram arising from $[\eta,[\eta,H]_{1}]_{2}$ which is
approximately equal to (b). In these diagrams, vertices of $H$ are represented
by circles, while vertices of $\eta$ are represented by squares.

In effect, the IMSRG(2∗) scheme approximates the goose-tank diagram in Fig.
3(b) with a second copy of the diagram in Fig. 3(d). This replacement is exact
for the fourth-order energy contribution, but only approximate for higher-
order terms. Explicitly, we have

$\displaystyle\Gamma^{(b)}_{ijkl}$
$\displaystyle\sim\tfrac{1}{2}\sum_{abcd}\bar{n}_{a}\bar{n}_{b}n_{c}n_{d}\frac{\Gamma_{ijdk}\Gamma_{cdab}\Gamma_{abcl}}{\Delta_{cdab}\Delta_{abcl}}$
(19) $\displaystyle\Gamma^{(d)}_{ijkl}$
$\displaystyle\sim\tfrac{1}{2}\sum_{abc}\bar{n}_{a}\bar{n}_{b}n_{c}n_{d}\frac{\Gamma_{ijdk}\Gamma_{cdab}\Gamma_{abcl}}{\Delta_{cdab}\Delta_{ijdk}}$

where $i,j$ are particle orbits and $k,l$ are holes. Diagram (d) comes with a
minus sign in (16) due to the commutator, so IMSRG(2∗) approximately captures
diagram (b) assuming $\Delta_{ijdk}\approx-\Delta_{cdab}$, which is true to
the extent that all the two-particle-two-hole denominators
$\Delta_{pp^{\prime}hh^{\prime}}$ are equal. Because these goose-tank
topologies are of pphh form, they feed back into the generator $\eta$ and have
a non-linear effect beyond the fourth-order energy. As we will see below, in
some difficult cases in which the IMSRG(2) flow diverges, the IMSRG(2*)
modification stabilizes the flow. An alternative strategy, in which the
relevant double commutators are evaluated exactly is presented in [40].

###  III.2 Perturbative triples

The remaining fourth order energy diagrams, the triples, are produced in the
IMSRG by terms involving the three-body generator [34, 17]. After solving the
flow equation (6) at the IMSRG(2) level, the transformed Hamiltonian is

$\bar{H}=e^{\Omega}He^{-\Omega}=\bar{H}_{\rm IMSRG(2)}+\bar{W}$ (20)

where $\bar{W}$ is the induced three-body interaction. We can directly
evaluate the contribution of this term to the ground state energy in second-
order perturbation theory. However, it is more convenient to use the Magnus
formulation here, because it makes generalizations beyond perturbative
counting and to the VS-IMSRG more straightforward. We wish to suppress the
off-diagonal part of $\bar{W}$, and the leading-order estimate of this, using
the White generator, is

$\bar{\Omega}=\frac{\bar{W}^{\rm od}}{\bar{\Delta}}.$ (21)

We perform a second unitary transformation on the Hamiltonian yielding
$\bar{\bar{H}}=e^{\bar{\Omega}}\bar{H}e^{-\bar{\Omega}}$. The zero-body piece
of this is

$\bar{\bar{H}}_{\rm 0b}=\bar{H}_{\rm 0b}+[\bar{\Omega},\bar{H}]_{\rm
0b}+\tfrac{1}{2}[\bar{\Omega},[\bar{\Omega},\bar{H}]]_{\rm 0b}+\ldots$ (22)

and the elipses indicate terms of order $\bar{\Omega}^{3}$ which we neglect
here. Assuming the denominator $\bar{\Delta}$ is chosen according the
partitioning we use for perturbative counting (i.e. Moller-Plesset or Epstein-
Nesbet), then taking just the zero-order piece of $\bar{H}$ gives
$[\bar{\Omega},\bar{H}_{\rm 1b,2b}^{[0]}]=-\bar{W}$ so

$\bar{\bar{H}}_{\rm 0b}=\bar{H}_{\rm
0b}+\tfrac{1}{2}[\bar{\Omega},\bar{W}]_{\rm 0b}+\ldots$ (23)

that is, we should take half the value of the $[3,3]_{0}$ commutator. Taking
the explicit form of the commutator, we have the perturbative correction to
the energy [34]

$\Delta
E_{[3]}=\frac{1}{(3!)^{2}}\sum_{abcijk}\bar{n}_{a}\bar{n}_{b}\bar{n}_{c}n_{i}n_{j}n_{k}\frac{\bar{W}_{ijkabc}\bar{W}_{abcijk}}{\bar{\Delta}_{ijkabc}}.$
(24)

There are several possible ways to approximate $\bar{W}$ and $\bar{\Delta}$
[34] which only involve a single $N^{7}$ commutator evaluation555More
precisely, this is $N_{h}^{3}N_{p}^{4}$, which is better than the worst
IMSRG(2) commutators., and which differ in the content beyond fourth order. A
strictly fourth-order correction would be
$\overline{W}=[\eta(0),\Gamma(0)]_{\rm 3b}$, with $\bar{\Delta}$ taken from
$H(0)$. In the Magnus formulation, we can just as easily take the full
$\Omega(\infty)$ obtained from the IMSRG(2) or IMSRG(2∗) solution. In fact, we
can also include higher nested commutators by using an intermediate

$\bar{W}=\sum_{k=1}^{\infty}\frac{1}{k!}[\Omega,H]^{(k)}_{\rm
3b}=[\Omega,\widetilde{H}]_{\rm 3b}$ (25)

with

$\widetilde{H}\equiv\sum_{k=0}^{\infty}\frac{1}{(k+1)!}[\Omega,H]^{(k)}.$ (26)

In this paper we will take the definition (25) and Møller-Plesset energy
denominators obtained from the one-body part of $H(\infty)$. We will indicate
this approximation Magnus(2∗)[3].

In the valence-space formulation, the connection to perturbation theory is
formally less straightforward, but (23) remains valid. We can still use (24)
to approximately capture the leading effects of the induced three-body
interaction. This effectively performs second-order perturbation theory on the
(ensemble) reference state. A more formally correct choice would be to use the
off-diagonal induced three-body interaction which connects valence
configurations to non-valence configurations, namely $W^{\rm
od}\in\\{W_{pppccc},W_{ppqccv},W_{ppqcvv},W_{ppqvvv}\\}$, where $c,v,q$
indicate core, valence, and excluded orbits, respectively, and
$p\in\\{q,v\\}$. The part of $\bar{\bar{H}}$ acting purely in the valence
space, including the induced 3N interaction, should then be included in the
valence space diagonalization (possibly perturbatively). The restriction of
external lines to valence orbits would keep the calculation tractable. This
direction will be left for a future study.

###  III.3 Truncation to $N^{7}$ scaling terms

Alternatively, we may assume that not all topologies contributing to the full
IMSRG(3) commutator expression are equally important. A pragmatic way to
approximate the full expression is to retain only those terms which scale as
$N^{7}$ or better. Explicitly, this amounts to neglecting the following
topologies shown in Fig. 1: $[3,3]_{2}$, $[3,2]_{3}$, and $[3,3]_{3}$. In this
approximation, errors made in the zero-body piece of $H(s)$ will be at least
fifth order in the potential, errors to the one-body and two-body parts will
be at least fourth order, and the error in the three-body part will be at
least third order. This truncation scheme was previously explored in Refs.
[41, 35], and we will refer to it as IMSRG(3N7). In principle we can also
consider terms scaling as $N^{8}$ or better, but in practice we find these
extra terms are not worth the substantial computational effort required for
direct evaluation.

##  IV Benchmarking

Before considering a more realistic case, we explore exactly-solvable
simplified models which separately probe the short-range and long-range
correlations present in a nucleus.

###  IV.1 Two fermions with a contact interaction

The simplest non-trivial many-body system consists of two particles. As
discussed in Section II.3, IMSRG(2) and even IMSRG(3) are not exact for this
system, because four-quasiparticle excitations are possible. Therefore, we can
gain some insight by considering a toy problem of two neutrons in a harmonic
trap with a contact interaction. We work in a basis spanned by three major
harmonic oscillator shells (i.e. $0s,0p,1s0d$). The Hamiltonian is

$H=\sum_{i}\epsilon_{i}a^{\dagger}_{i}a_{i}+\tfrac{1}{4}\sum_{ijkl}V_{ijkl}a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}$
(27)

where the single-particle energy is $\epsilon_{i}=2n_{i}+\ell$ (we measure
energy in terms of the oscillator energy $\hbar\omega$) and the potential is

$V(\vec{r})=\frac{g}{({2\pi})^{3/2}}\delta(\vec{r}).$ (28)

The relative coordinate $\vec{r}$ is expressed in units of the oscillator
length $b=\sqrt{\hbar/m\omega}$. With this definition, $g<0$ corresponds to an
attractive interaction, and the normalization is chosen so that $\langle
0s0s|V|0s0s\rangle=g$.

While we must retain induced three-body (and four-body) terms in $H$ to obtain
the exact result, there are no three-particle or three-hole configurations
possible for an $A=2$ system. Only matrix elements involving $pph$ and $hhp$
configurations need to be considered. Three-body terms in the generator, which
are of the form $\eta_{ppphhh}$, will not contribute. Consequently, the
perturbative triples discussed in section III.2 give no contribution to the
energy, and the IMSRG(2*) approximation is exact through $g^{4}$ in
perturbation theory.

Figure 4: Correlation energy of two identical fermions in a harmonic trap with
a contact interaction, obtained in various approximation schemes as a function
of the coupling strength $g$. The lower panel shows the error in the energy as
a percent of the correlation energy.

Fig. 4 displays the ground state energy as a function of $g$, computed at
various levels of approximation. We include perturbation to second and third
order, labeled MBPT(2) and MBPT(3), respectively, as well as MBPT(3) plus
particle-particle and hole-hole ladders resummed to all orders. We also
include IMSRG(2), IMSRG(2*), IMSRG(3N7), and the full IMSRG(3). In the top
panel of Fig. 4, we observe that for large positive coupling $g$, perturbation
theory breaks down completely and resummed ladders do not fix the problem. For
attractive couplings, the relevant physics is captured by the ladder
approximation. In contrast, the IMSRG(2) is accurate for repulsive couplings,
but diverges for sufficiently attractive coupling. All approximations beyond
IMSRG(2) are sufficiently accurate that they cannot be distinguished in the
top panel. In the bottom panel of Fig. 4, we show the error relative to the
exact solution, zoomed in to visualize the difference between the IMSRG
approximations. It appears that there is little further improvement in going
from IMSRG(2*) to full IMSRG(3).

Figure 5: Flow of (a) zero-body term $E_{0}$ and (b) off-diagonal part of the
two-body vertex $\Gamma_{pphh}$, for the two-fermion problem with $g=-2$,
computed with the IMSRG(2) and IMSRG(2∗) approximations.

For $g\lesssim-1.8$, the IMSRG(2) flow diverges, while the flow converges to
the fixed point for the IMSRG(2∗) and other approximations of IMSRG(3), as
illustrated in Fig. 5. Fig. 5(b), shows the flow of the off-diagonal part of
the two-body vertex $\Gamma_{pphh}$. With the IMSRG(2) approximation, the off-
diagonal pieces are not suppressed, meaning $\frac{d}{ds}|\Gamma_{pphh}|>0$.
This can be traced back to the contribution illustrated in Fig. 6 in which the
two-body vertex in the particle-particle channel
$\Gamma_{ppp^{\prime}p^{\prime}}$ contracts with the generator
$\eta_{p^{\prime}p^{\prime}hh}$ to give a contribution to
$\frac{d}{ds}\Gamma_{pphh}$. The contact interaction is approximately
separable, i.e. $\Gamma_{abcd}\sim gF_{ab}F_{cd}$ for some one-body matrix
$F$. This leads to a coherent effect

$\displaystyle\frac{d}{ds}\left(gF_{pp}F_{hh}\right)$
$\displaystyle\sim-\sum_{p^{\prime}p^{\prime}}\left(gF_{pp}F_{p^{\prime}p^{\prime}}\right)\left(gF_{p^{\prime}p^{\prime}}F_{hh}\right)$
(29)
$\displaystyle\sim(gF_{pp}F_{hh})\cdot(-g)\cdot\sum_{p^{\prime}p^{\prime}}(F_{p^{\prime}p^{\prime}})^{2}.$

For $g<0$ the contribution to $\frac{d}{ds}\Gamma_{pphh}$ has the same sign as
$\Gamma_{pphh}(s)$, driving an exponential enhancement.

Figure 6: Diagram illustrating the commutator term
$[\eta,\Gamma_{ppp^{\prime}p^{\prime}}]$ which contributes to the increasing
magnitude of the off-diagonal two-body vertex $\Gamma_{pphh}$.

In contrast, the diagrams in Fig. 3 give a contribution with the opposite sign
of $\Gamma_{pphh}$, leading to suppression. This can be seen by considering
the expression for the intermediate $\chi_{hh}$ as defined in (16). For the
White generator this becomes

$\displaystyle\chi_{hh}$
$\displaystyle=\sum_{p^{\prime}}(\eta_{hhp^{\prime}p^{\prime}}\Gamma_{p^{\prime}p^{\prime}hh}-\Gamma_{hhp^{\prime}p^{\prime}}\eta_{p^{\prime}p^{\prime}hh})$
(30)
$\displaystyle=-2\sum_{p^{\prime}}\frac{|\Gamma_{p^{\prime}p^{\prime}hh}|^{2}}{\Delta_{p^{\prime}p^{\prime}hh}}$

which is manifestly negative. Contracting $\chi_{hh}$ with $\eta_{pphh}$ then
yields a contribution $\frac{d}{ds}\Gamma_{pphh}\sim-\Gamma_{pphh}$, leading
to suppression.

This illustrates how, in addition to restoring the fourth-order quadruples
contribution to the ground state energy, the IMSRG(2∗) approximation—and the
more expensive IMSRG(3N7) and IMSRG(3)—stabilizes the IMSRG flow in the non-
perturbative regime. For this toy problem, this constitutes the main effect of
IMSRG(3).

Figure 7: Error as a function of coupling constant $g$ for 8 identical
fermions in a harmonic trap.

In order to illustrate the impact of the perturbative triples correction, we
repeat the calculation with 8 particles in the trap. The error as a function
of the coupling $g$ is shown in Fig. 7. The behavior is qualitatively similar
to the two-particle case, with the exception that the perturbative triples now
have an effect. For repulsive $g>0$, the perturbative triples correction is
substantial and accounts for nearly all of the difference between IMSRG(2∗)
and IMSRG(3N7). On the other hand, for $g<0$, IMSRG(2∗) overbinds relative to
IMSRG(3N7); it is clear that the perturbative triples correction, which
necessarily lowers the energy, cannot correct for this difference.

###  IV.2 Lipkin-Meshkov-Glick model

Another toy problem which can illuminate many-body approximation schemes is
the Lipkin-Meshkov-Glick (LGM) model [42]. The model consists of $N$ particles
distributed across two $N$-fold degenerate levels separated by energy
$\epsilon$. The Hamiltonian is

$H=\tfrac{1}{2}\epsilon\sum_{p\sigma}\sigma
a^{\dagger}_{p\sigma}a_{p\sigma}+\tfrac{1}{2}V\sum_{pp^{\prime}\sigma}a^{\dagger}_{p\sigma}a^{\dagger}_{p^{\prime}\sigma}a_{p^{\prime}-\sigma}a_{p-\sigma},$
(31)

where $\sigma\in\\{-1,1\\}$ labels the lower and upper levels, and
$p,p^{\prime}$ run from 1 to $N$, labeling the different degenerate states
within a level.

This model can be solved analytically with a quasi-spin formulation,
introducing the operators

$K_{\pm}=\sum_{p}a^{\dagger}_{p\pm 1}a_{p\mp 1}\hskip 10.00002pt,\hskip
10.00002ptK_{z}=\tfrac{1}{2}\sum_{p\sigma}\sigma
a^{\dagger}_{p\sigma}a_{p\sigma}$ (32)

which obey the usual angular momentum commutation rules (see Appendix C). In
terms of the quasi-spin operators, the Hamiltonian (31) can be expressed as

$H=\epsilon K_{z}+\tfrac{1}{2}V(K_{+}^{2}+K_{-}^{2}).$ (33)

$H$ commutes with the operator
$K^{2}=K_{z}^{2}+\tfrac{1}{2}(K_{+}K_{-}+K_{-}K_{+})$, which has eigenvalues
$k(k+1)$, with $k\leq\frac{N}{2}$. In the thermodynamic ($N\to\infty$) limit,
the LGM model exhibits a phase transition from an ungapped ferromagnetic phase
for $NV/\epsilon<1$ to a gapped unmagnetized phase for $NV/\epsilon>1$.

We focus here on the specific case $N=8$, which has the analytical
solution666There appears to be a typo in the solution presented in Ref. [42].
for the ground state energy [42]

$\frac{E}{\epsilon}=-\sqrt{10+118\left(\frac{V}{\epsilon}\right)^{2}+6\sqrt{1-2\left(\frac{V}{\epsilon}\right)^{2}+225\left(\frac{V}{\epsilon}\right)^{4}}}.$
(34)

The LGM model has been solved in coupled cluster theory, where it has an
analytical result in the doubles (CCD) approximation [43] (see also [44]). The
CCD correlation energy is

$\frac{E_{\rm
corr}}{\epsilon}=-\frac{N(N-1)}{2\alpha}\left(1-\sqrt{1-\alpha\frac{V^{2}}{\epsilon^{2}}}\right)$
(35)

where $\alpha\equiv N^{2}-7N+9$. We can see that for $N\geq 6$, $\alpha>0$ and
the energy becomes complex for sufficiently large values of $V/\epsilon$. For
the specific case $N=8$ which we consider here, CCD fails for
$NV/\epsilon>\sqrt{64/17}\approx 1.94$. For comparison with the IMSRG(2) we
also consider an approximation of CCD in which we omit the intermediate three-
body operator arising in the CCD decoupling condition [37]. We obtain the same
expression for the correlation energy with the modification
$\alpha\to\tilde{\alpha}\equiv N^{2}-5N+7$. Note that because this model
exhibits an “excitation-parity” symmetry, single and triple excitations have
no effect, so that CCD is equivalent to CCSDT.

Moving on to the IMSRG solution, we note that the LGM model has been
frequently studied within flow equation approach [45, 46, 31, 47, 48, 49],
notably by Pirner and Friman [45], who used the quasi-spin operators (without
normal ordering) to formulate the flow equation. They found it necessary to
include an induced three-body term, but even then obtained results which were
less accurate than second order perturbation theory. More recently, the flow
equation approach was applied to obtain $1/N$ corrections to the thermodynamic
limit [49]. Here we are interested in finite system, and in understanding the
impact of truncations in the particle-rank of the flowing operators. Making an
approximation equivalent to IMSRG(2) requires the use of normal-ordered
operators. Commutators of normal-ordered strings of quasi-spin operators are
presented in Appendix C. While we may construct the IMSRG flow equations in
terms of the normal-ordered quasi-spin operators, integrating these equations
must still be done numerically. However, in the Magnus formulation, we may
directly solve for the Magnus operator at the fixed point, obtaining (see
Appendix C)

$\tanh\left(\sqrt{4\Omega_{\rm
2b}^{2}\tilde{\alpha}}\right)=\sqrt{\tilde{\alpha}}\frac{V}{\epsilon}$ (36)

and an energy which is identical to the CCD energy with
$\alpha\to\tilde{\alpha}$. Clearly the IMSRG(2) will fail if
$\sqrt{\tilde{\alpha}}V/\epsilon>1$.

We now turn to the effects of IMSRG(3). Due to the excitation-parity symmetry
of the LGM model, the three-body pieces of $\eta$ and $\Omega$ vanish.
Consequently, the perturbative triples correction described in section III.2
has no effect. The only way in which induced three-body operators affect the
ground state energy by modifying the flow for $V(s)$ which then subsequently
gets contracted with $\Omega$ to generate a contribution to the zero-body
term.

The factor of $\tilde{\alpha}$ arises because the double commutator
$[K_{\pm}^{2},[K_{+}^{2},K_{-}^{2}]]=\pm 8\tilde{\alpha}K_{\pm}^{2}$ in the
NO2B approximation. If we also include in our set of flowing operators a term
proportional to the three-body operator $\\{K_{+}K_{z}K_{-}\\}$, we capture
the goose-tank diagrams of Fig. 3, and the double-commutator is modified so
that $\tilde{\alpha}\to\alpha$, exactly reproducing the CCD energy. In the
IMSRG(2∗) approximation, we modify the double commutator by adding the
commutator $[\Omega,\chi]$ with the auxiliary operator $\chi$, which also has
the effect of replacing $\tilde{\alpha}\to\alpha$, without the need to
explicitly construct a three-body operator.

There are two additional Hermitian three-body operators consistent with the
symmetries of the problem: $\\{K_{z}^{3}\\}$ and
$(\\{K_{+}^{2}K_{z}\\}+\\{K_{z}K_{-}^{2}\\})$. Including these spoils the
simple structure obtained in the IMSRG(2) approximation, so an analytic
solution is less straightforward. However, it is easy to include them
numerically. The results are shown in Fig. 8.

Figure 8: Energy for the $N=8$ LGM model with various approximate methods. See
text for details.

As mentioned above, the IMSRG(2*) energy is identical to the CCD energy, while
IMSRG(2) is identical to CCD with the intermediate three-body operator
removed, indicated CCD′ in Fig. 8. IMSRG(3N7) does better than both of these,
and IMSRG(3) gives a significant further improvement. As mentioned above, due
to the excitation-parity symmetry of the problem, there is no improvement from
CCSD to CCSDT. Even in the IMSRG(3) approximation the generator $\eta$ and the
Magnus operator $\Omega$ are purely two-body, which means that, as in CCD,
both the IMSRG(2) and IMSRG(3) transformations are specified by a single
parameter, $\Omega_{\rm 2b}$. All the improvement from IMSRG(2) to IMSRG(3)
comes from more accurately evaluating the transformation (8) by including
intermediate three-body operators which subsequently get contracted back down
to lower-body operators. In the specific case in which the three-body operator
is immediately contracted back down, e.g. $[\Omega,[\Omega,H]_{\rm 3b}]_{\rm
2b}$, it is possible to factorize the expression for the double nested
commutator so that it scales as $N^{6}$ [50, 34, 7]. For the LGM model,
including this factorizable correction for all nested commutators in (8)
exactly reproduces the IMSRG(3N7) results. An application to the more
realistic case will be presented in a separate work [40].

Thanks to the simplicity of the LGM model, we may also evaluate the
transformation (8) directly in the $A$-body space using the quasispin basis.
For $N=8$ this involves the manipulation of $5\times 5$ matrices. Taking the
ansatz $\Omega=\tfrac{1}{2}\Omega_{2b}(K_{+}^{2}-K_{-}^{2})$, we solve the
flow equations numerically without truncating any intermediate operators. The
result is labeled “$\Omega_{2b}$ Exact” in Fig. 8. It lies above the full
exact result, as it satisfies the variational principle.

Finally, because the flowing Hamiltonian has essentially two parameters in the
IMSRG(2) and IMSRG(2∗) approximations, we can visualize the SRG trajectory, as
shown in Fig. 9. (The symmetry about the $x$ axis reflects the symmetry of the
LGM model under $V\to-V$.) Well below the critical coupling strength
$\sqrt{\tilde{\alpha}}V/\epsilon=1$, the IMSRG(2) and IMSRG(2∗) are both
accurate and follow essentially the same trajectory. Near the critical
coupling, the IMSRG(2) diverges significantly from the IMSRG(2∗) trajectory,
and beyond the critical coupling the fixed point at $V=0$ disappears, so the
flow continues to the fixed point at infinity.

Figure 9: IMSRG flow in the LGM model, for the IMSRG(2) and IMSRG(2∗)
approximations for several values of $V(0)$. Arrows indicate the direction of
the flow. The red lines indicate the exact solutions from (34). cf Ref. [48].

To be clear, the observations of this section do not mean that the three-body
part of the Magnus operator is unimportant in realistic applications; it is
identically zero here because of a symmetry of the LGM model which is not
present in general777Though, for closed harmonic oscillator shells, the usual
parity symmetry does impart an approximate excitation-parity symmetry to low-
lying states.. However, it does demonstrate the significant improvements that
can be obtained in the IMSRG by restoring many-body terms inside nested
commutators.

###  IV.3 Realistic nuclear interaction with the valence-space IMSRG

Figure 10: Ground state energy and excitation energy of the first $2^{+}$
state in carbon isotopes, computed in an oscillator space with $e_{max}=3$,
$\hbar\omega=16$ MeV.

We next turn to a more realistic nuclear structure calculation, in which we
consider the even-$A$ carbon isotopes with $A=8-16$ computed with a pure two-
body interaction derived from chiral EFT with a cutoff $\Lambda=500$ MeV and
no free-space SRG softening. We use this relatively hard interaction to
emphasize non-perturbative effects. Our current implementation of the IMSRG(3)
limits us to $e_{\rm max}\lesssim 4$, so we cannot perform fully-converged
calculations to compare with experiment. (In any case, comparing with
experiment would conflate many-body errors with errors due to the input
nuclear force.) Instead, we compare with truncated configuration interaction
(CI) in a given $e_{\rm max}$ space. This introduces a contamination due to
spurious excitations of the center of mass [51, 27], but we find that for this
specific system the contamination is negligible by $e_{\rm max}=3$.

Ground state energies and $2^{+}_{1}$ excitation energies of some even-$A$
carbon isotopes are shown in Fig. 10 for $e_{\rm max}=3$. To focus on the
effect of the IMSRG evolution, for the ground state energies we subtract off
the energy obtained by diagonalizing in the valence space with the unevolved
$s=0$ Hamiltonian. Various levels of approximation, from no IMSRG evolution
(labeled “VS-HF”, since we use a Hartree-Fock basis), to full IMSRG(3), are
compared with the truncated configuration-interaction (CI) calculations
extrapolated to the full CI (FCI) result. We truncate on $N_{\rm max}$, i.e.
the maximum number of harmonic oscillator quanta of excitation, up to $N_{\rm
max}=$ 6 or 8. The gray error bands in Fig. 10 indicate the difference from
the largest $N_{\rm max}$ calculation performed to the $N_{\rm max}\to\infty$
extrapolated value. The first observation is that going from IMSRG(2) to
IMSRG(3) generally improves agreement with FCI for both the ground state
energy and the excitation energy. We also notice that the IMSRG(3N7)
approximation accurately reproduces the full IMSRG(3) result. The IMSRG(2∗)
approximation overbinds in the valence space formulation and slightly
underbinds in the single-reference formulation. Nevertheless, it is notable
that the SR-IMSRG(2∗) flow converges for 12C, while it diverges for SR-
IMSRG(2).

The IMSRG(2) approximation consistently yields excitation energies that are
too high, while this is corrected by going to IMSRG(3N7). It appears that the
most important modification to the effective valence space interaction in
going from IMSRG(2) to IMSRG(3N7) is that the off-diagonal matrix element
$\langle p_{3/2}p_{3/2}|V|p_{1/2}p_{1/2}\rangle_{J=0}$ is reduced in magnitude
by $\sim\\!0.5$ MeV. This term, which is too large in the IMSRG(2)
approximation, drives 2p-2h correlations which provide too much correlation
energy for the ground state. These correlations are less active in the $2^{+}$
excited state, and so the net result is an over-predicted excitation energy.
It is also interesting to note that the IMSRG(2*) approximation has
essentially no impact on these off-diagonal valence space matrix elements, and
thus gives no improvement for the $2^{+}$ excitation energy. This can be
understood because the IMSRG(2*) approximation essentially generates the
goose-tank diagrams of Fig. 3 via a contraction of a one-body operator $\chi$
with the two-body part of the generator, so that the main effect is to modify
the off-diagonal part of $H$. Valence-to-valence matrix elements are not
included in the definition of off-diagonal and so this term receives
relatively little modification.

##  V Computational cost

Finally, we briefly discuss the computational cost of our current
implemententation of the IMSRG(3), and approximations. While this is to some
extent specific to this implementation, other implementations will share the
general features, and this also helps to illuminate where future improvements
should be focused. All commutator expressions involving a three-body operator
have been recast as matrix matrix multiplications, leading to significant
speedup with the exception of the $[2,2]\to 3$ expression. In this case, we
decided that the modest speedup was not worth the considerable complication in
the code. Multi-threaded parallelism via OpenMP is utilized throughout these
routines. The time for each expression is shown in Fig. 11. Which commutator
expression is the most expensive depends on the model space and system
studied, but beyond the smallest spaces the $[2,2]\to 3$ and $[2,3]\to 2$
expressions are the main bottleneck of IMSRG(3N7) calculations. For full
IMSRG(3), terms with intermediate pph lines dominate, due to the angular
momentum recoupling required to cast them as matrix multiplications. We
observe that the naive scaling estimate is indeed naive; the $[3,3]\to 3$
commutator with ppp or hhh intermediate states scales as $N^{9}$, but it takes
a negligible amount of time because it is highly amenable to matrix
multiplication without any required recouplings or antisymmetrizations.

Figure 11:  Compute time per commutator on a single node with 64 threads, for
three levels of approximation. The labels indicate the particle rank of the
two operators entering the expression and the resulting operator. The
following strings of p’s and h’s indicate whether the contracted indices are
particle or hole, for cases where the routine was split up according to this
feature.

##  VI Conclusion

We have performed a detailed exploration of the impact of retaining three-body
operators within the IMSRG, building on previous work [35]. While a full
inclusion of three-body operators remains too expensive for practical
calculations, we find that more manageable approximations capture the most
important corrections. We explored two toy models, selected to emphasize
short-range and long-range correlations, respectively. We find that in the
case of a short-range interaction, the main error made in the IMSRG(2) is
corrected by including the effects of the goose-tank diagrams in Fig. 3, which
can be done without changing the computational scaling. We also investigated
the origin of the diverging IMSRG flow for sufficiently large coupling
constant $g$, and related it to coherent effects due to the approximately
separable nature of the interaction; these coherent effects were offset by
including the missing goose-tank diagrams, stabilizing the flow.

In the case of long-range correlations, as captured by the LGM model, we
presented an analytic solution for the IMSRG(2) approximation, and showed that
including intermediate three-body can lead to a significant improvement in
accuracy over coupled cluster. Finally, we have explored the impact of three-
body operators in the valence-space IMSRG with a realistic interaction. We
found that the systematic over-prediction of $2^{+}$ excitation energies in
the IMSRG(2) approximation is essentially corrected by including three-body
operators during the IMSRG evolution, and that this improvement is achieved
already at the IMSRG(3N7) level of approximation. Finally, we have presented
arguments, supported by our calculations, that the most important effect of
three-body operators arises as intermediate terms in nested commutators. The
leading corrections of this type can be factorized to scale as $N^{6}$, and so
may be employed in large-scale calculations. This is presented in a separate
paper [40].

###### Acknowledgements.

S.R.S. would like to thank Matthias Heinz for help on optimizing the
implementation and finding a typo in the commutator expressions. S.R.S. and
B.C.H. acknowledge support from the U.S. Department of Energy (DoE), Office of
Science, under SciDAC-5 (NUCLEI collaboration) under grant DE-FG02-97ER41014.
This research was supported in part by the Notre Dame’s Center for Research
Computing. This work benefitted from the use of the amc code [52] to validate
the $J$-coupled expressions. The imsrg++ code [53] employed in this work uses
the Armadillo library [54, 55].

##  Appendix A Uncoupled commutator expressions

We need the expression for the commutator $C=[A,B]$ where $A$, $B$, and $C$
are all Fock-space operators including up to 3 body operators, e.g.

$\displaystyle C=C_{0}$
$\displaystyle+\sum_{ij}C_{ij}\\{a^{\dagger}_{i}a_{j}\\}+\frac{1}{4}\sum_{ijkl}C_{ijkl}\\{a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}\\}$
(37)
$\displaystyle+\frac{1}{36}\sum_{ijklmn}C_{ijklmn}\\{a^{\dagger}_{i}a^{\dagger}_{j}a^{\dagger}_{k}a_{n}a_{m}a_{l}\\}$

and likewise for $A$ and $B$. In the following, the permutation operator
$P_{ij}$ exchanges the indices $i$ and $j$ of the expression to its right,
with no implicit minus sign. For example
$P_{ij}A_{bjal}B_{aibk}=A_{bial}B_{ajbk}$. We also use the operator

$P_{ij/k}\equiv 1-P_{ik}-P_{jk}.$ (38) $\displaystyle C_{0}=$
$\displaystyle\sum_{ab}(n_{a}-n_{b})A_{ab}B_{ba}+\frac{1}{4}\sum_{abcd}n_{a}n_{b}\bar{n}_{c}\bar{n}_{d}(A_{abcd}B_{cdab}-B_{abcd}A_{cdab})$
(39)
$\displaystyle+\frac{1}{36}\sum_{abcdef}n_{a}n_{b}n_{c}\bar{n}_{d}\bar{n}_{e}\bar{n}_{f}(A_{abcdef}B_{defabc}-B_{abcdef}A_{defabc})$
$\displaystyle C_{ij}=$
$\displaystyle\sum_{a}(A_{ia}B_{aj}-B_{ia}A_{aj})+\sum_{ab}(n_{a}-n_{b})(A_{ab}B_{biaj}-B_{ab}A_{biaj})$
(40)
$\displaystyle+\frac{1}{2}\sum_{abc}(n_{a}n_{b}\bar{n}_{c}+\bar{n}_{a}\bar{n}_{b}n_{c})(A_{ciab}B_{abcj}-B_{ciab}A_{abcj})$
$\displaystyle+\frac{1}{4}\sum_{abcd}n_{a}n_{b}\bar{n}_{c}\bar{n}_{d}(A_{abcd}B_{cdiabj}+A_{abicdj}B_{cdab}-B_{abcd}A_{cdiabj}-B_{abicdj}A_{cdab})$
$\displaystyle+\frac{1}{12}\sum_{abcde}(n_{a}n_{b}\bar{n}_{c}\bar{n}_{d}\bar{n}_{e}+\bar{n}_{a}\bar{n}_{b}n_{c}n_{d}n_{e})(A_{abicde}B_{cdeabj}-B_{abicde}A_{cdeabj})$
$\displaystyle C_{ijkl}=$
$\displaystyle\sum_{a}(A_{ia}B_{ajkl}+A_{ja}B_{iakl}-B_{ijal}A_{ik}-B_{ijka}A_{al})$
(41)
$\displaystyle-\sum_{a}(B_{ia}A_{ajkl}+B_{ja}A_{iakl}-A_{ijal}B_{ka}-A_{ijka}B_{ak})$
$\displaystyle+\frac{1}{2}\sum_{ab}(\bar{n}_{a}\bar{n}_{b}-n_{a}n_{b})(A_{ijab}B_{abkl}-B_{ijab}A_{abkl})$
$\displaystyle-\sum_{ab}(n_{a}-n_{b})(1-P_{ij})(1-P_{kl})A_{bjal}B_{aibk}$
$\displaystyle+\sum_{ab}(n_{a}-n_{b})(A_{ab}B_{ijbkla}-B_{ab}A_{ijbkla})$
$\displaystyle+\frac{1}{2}\sum_{abc}(n_{a}n_{b}\bar{n}_{c}+\bar{n}_{a}\bar{n}_{b}n_{c})(1-P_{ij})(A_{icab}B_{ajbklc}-B_{icab}A_{abjklc})$
$\displaystyle-\frac{1}{2}\sum_{abc}(n_{a}n_{b}\bar{n}_{c}+\bar{n}_{a}\bar{n}_{b}n_{c})(1-P_{kl})(B_{ijcabl}A_{abkc}-A_{ijcabl}B_{abkc})$
$\displaystyle+\frac{1}{6}\sum_{abcd}(n_{a}n_{b}n_{c}\bar{n}_{d}-\bar{n}_{a}\bar{n}_{b}\bar{n}_{c}n_{d})(A_{ijdabc}B_{abckld}-B_{ijdabc}A_{abckld})$
$\displaystyle+\frac{1}{4}\sum_{abcd}(\bar{n}_{a}\bar{n}_{b}n_{c}n_{d}-n_{a}n_{b}\bar{n}_{c}\bar{n}_{d})(1-P_{ij})(1-P_{kl})A_{abicdk}B_{cdjabl}$
$\displaystyle C_{ijklmn}=$
$\displaystyle\sum_{a}P_{i/jk}(A_{ia}B_{ajklmn}-B_{ia}A_{ajklmn})-P_{l/mn}(B_{ijkamn}A_{al}-A_{ijkamn}B_{al})$
(42)
$\displaystyle-\sum_{a}P_{ij/k}P_{lm/n}(A_{ijna}B_{kalm}-B_{ijna}A_{kalm})$
$\displaystyle+\frac{1}{2}\sum_{ab}(\bar{n}_{a}\bar{n}_{b}-n_{a}n_{b})P_{ij/k}(A_{ijab}B_{abklmn}-B_{ijab}A_{abklmn})$
$\displaystyle-\frac{1}{2}\sum_{ab}(\bar{n}_{a}\bar{n}_{b}-n_{a}n_{b})P_{lm/n}(B_{ijkabn}A_{ablm}-A_{ijkabn}B_{ablm})$
$\displaystyle-\sum_{ab}(n_{a}\bar{n}_{b}-\bar{n}_{a}n_{b})P_{ij/k}P_{lm/n}(A_{bkan}B_{ijalmb}-B_{bkan}A_{ijalmb})$
$\displaystyle+\frac{1}{6}\sum_{abc}(n_{a}n_{b}n_{c}+\bar{n}_{a}\bar{n}_{b}\bar{n}_{c})(A_{ijkabc}B_{abclmn}-B_{ijkabc}A_{abclmn})$
$\displaystyle-\frac{1}{2}\sum_{abc}(n_{a}n_{b}\bar{n}_{c}+\bar{n}_{a}\bar{n}_{b}n_{c})P_{ij/k}P_{lm/n}(A_{ijcabn}B_{abklmc}-B_{ijcabn}A_{abklmc})$

##  Appendix B Commutator expressions in $J$-coupled form

To take advantage of rotational symmetry, we work with $J$-coupled matrix
elements. To do this, we switch our notation so that we now explicitly
indicate the projection quantum number $m$ in the uncoupled basis. One body
terms are written as

$C_{im_{i}jm_{j}}=C_{ij}\delta_{m_{i}m_{j}}$ (43)

while the two and three-body terms are written in terms of unnormalized
$J$-coupled matrix elements

$C_{im_{i}jm_{j}km_{k}lm_{l}}=\sum_{JM}\mathcal{C}_{j_{i}m_{i}j_{j}m_{j}}^{JM}\mathcal{C}_{j_{k}m_{k}j_{l}m_{l}}^{JM}C^{J}_{ijkl}$
(44) $\displaystyle
C_{im_{i}jm_{j}km_{k}lm_{l}mm_{m}nm_{n}}=\sum_{\begin{subarray}{c}J_{1}J_{2}\mathcal{J}\\\
M_{1}M_{1}\mathcal{M}\end{subarray}}\mathcal{C}_{j_{i}m_{i}j_{j}m_{j}}^{J_{1}M_{1}}\mathcal{C}_{J_{1}M_{1}j_{k}m_{k}}^{\mathcal{J}\mathcal{M}}$
(45)
$\displaystyle\times\mathcal{C}_{j_{l}m_{l}j_{m}m_{m}}^{J_{2}M_{2}}\mathcal{C}_{J_{2}M_{2}j_{n}m_{n}}^{\mathcal{J}\mathcal{M}}C^{J_{1}J_{2}\mathcal{J}}_{ijklmn}.$

The $\mathcal{C}$ in (44) and (45) are Clebsch-Gordan coefficients. By
inserting (43) (44) and (45) into (39), (40), (41), and (42), we obtain the
$J$-coupled commutator expressions.

In order to make the expressions somewhat more compact, in the following we
use orbit labels to indicate the associated angular momentum, e.g. in 6$j$,
9$j$ symbols and phase factors so that $(-1)^{i+j-J}$ is shorthand for
$(-1)^{j_{j}+j_{j}-J}$. We also introduce the coupled permutation operators

$P^{J}_{ij}\equiv(-1)^{i+j-J}P_{ij}$ (46) $\displaystyle
P^{J_{1}\mathcal{J}}_{ij/k}\equiv 1$
$\displaystyle-\sum_{J_{1}^{\prime}}\langle
ijkJ_{1}\mathcal{J}|kjiJ_{1}^{\prime}\mathcal{J}\rangle
P_{ik}^{J_{1}\rightarrow J_{1}^{\prime}}$ (47)
$\displaystyle-\sum_{J_{1}^{\prime}}\langle
ijkJ_{1}\mathcal{J}|ikjJ_{1}^{\prime}\mathcal{J}\rangle
P_{jk}^{J_{1}\rightarrow J_{1}^{\prime}}$

The brackets in (47) are recoupling coefficients which may be expressed in
terms of 6$j$ symbols

$\langle
ijkJ_{1}\mathcal{J}|kjiJ_{1}^{\prime}\mathcal{J}\rangle=\hat{J}_{1}\hat{J}_{1}^{\prime}(-1)^{2i+2j+2k}\begin{Bmatrix}i&j&J_{1}\\\
k&\mathcal{J}&J_{1}^{\prime}\end{Bmatrix}$ (48)

and

$\langle
ijkJ_{1}\mathcal{J}|ikjJ_{1}^{\prime}\mathcal{J}\rangle=\hat{J}_{1}\hat{J}_{1}^{\prime}(-1)^{j+k+J_{1}+J_{1}^{\prime}}\begin{Bmatrix}j&i&J_{1}\\\
k&\mathcal{J}&J_{1}^{\prime}\end{Bmatrix}.$ (49)

The notation $J_{1}\rightarrow J_{1}^{\prime}$ in the superscript in (47)
indicates that every instance of $J_{1}$ to the right of the permutation
operator should be replaced with $J_{1}^{\prime}$. We use the usual notation
$\hat{J}\equiv\sqrt{2J+1}$. With this notation, we have for our
antisymmetrized matrix elements

$P^{J}_{ij}C^{J}_{ijkl}=-C^{J}_{ijkl}$ (50)

and

$P^{J_{1}\mathcal{J}}_{ij/k}C^{J_{1}J_{2}\mathcal{J}}_{ijklmn}=3~{}C^{J_{1}J_{2}\mathcal{J}}_{ijklmn}.$
(51)

As in previous presentations, we employ Pandya-transformed operators, denoted
with an overbar

$\bar{C}^{J}_{i\bar{j}k\bar{l}}\equiv-\sum_{J^{\prime}}\hat{J^{\prime}}^{2}\begin{Bmatrix}i&j&J\\\
k&l&J^{\prime}\end{Bmatrix}C^{J^{\prime}}_{ilkj}.$ (52)

We also define a transformed three-body operator as

$\displaystyle\bar{C}^{J^{\prime}}_{i\bar{l};(cdJ_{cd})(\overline{abJ_{ab}})}\equiv-\sum_{\mathcal{J}}$
$\displaystyle\hat{\mathcal{J}}^{2}(-1)^{i+\mathcal{J}}\begin{Bmatrix}i&l&J^{\prime}\\\
J_{cd}&J_{ab}&\mathcal{J}\end{Bmatrix}$ (53) $\displaystyle\times
C^{J_{ab}J_{cd}\mathcal{J}}_{abicdl}.$

With these definitions, the $J$-coupled commutator expressions become

$\displaystyle C_{0}=$
$\displaystyle\sum_{ab}\hat{j}^{2}_{a}(n_{a}-n_{b})A_{ab}B_{ba}+\frac{1}{4}\sum_{abcd}\sum_{J}\hat{J}^{2}n_{a}n_{b}\bar{n}_{c}\bar{n}_{d}(A^{J}_{abcd}B^{J}_{cdab}-B^{J}_{abcd}A^{J}_{cdab})$
(54)
$\displaystyle+\frac{1}{36}\sum_{abcdef}\sum_{J_{1}J_{2}\mathcal{J}}n_{a}n_{b}n_{c}\bar{n}_{d}\bar{n}_{e}\bar{n}_{f}\hat{\mathcal{J}}^{2}(A^{J_{1}J_{2}\mathcal{J}}_{abcdef}B^{J_{2}J_{1}\mathcal{J}}_{defabc}-B^{J_{1}J_{2}\mathcal{J}}_{abcdef}A^{J_{2}J_{1}\mathcal{J}}_{defabc})$
$\displaystyle C_{ij}=$
$\displaystyle\sum_{a}(A_{ia}B_{aj}-B_{ia}A_{aj})+\sum_{ab}\sum_{J}(n_{a}-n_{b})\frac{\hat{J}^{2}}{\hat{j}_{i}^{2}}(A_{ab}B^{J}_{biaj}-B_{ab}A^{J}_{biaj})$
(55)
$\displaystyle+\frac{1}{2}\sum_{abc}\sum_{J}(n_{a}n_{b}\bar{n}_{c}+\bar{n}_{a}\bar{n}_{b}n_{c})\frac{\hat{J}^{2}}{\hat{j}_{i}^{2}}(A^{J}_{ciab}B^{J}_{abcj}-B^{J}_{ciab}A^{J}_{abcj})$
$\displaystyle+\frac{1}{4}\sum_{abcd}\sum_{J\mathcal{J}}n_{a}n_{b}\bar{n}_{c}\bar{n}_{d}\frac{\hat{\mathcal{J}}^{2}}{\hat{j}_{i}^{2}}(A^{J}_{abcd}B^{JJ\mathcal{J}}_{cdiabj}+A^{JJ\mathcal{J}}_{abicdj}B^{J}_{cdab}-B^{J}_{abcd}A^{JJ\mathcal{J}}_{cdiabj}-B^{JJ\mathcal{J}}_{abicdj}A^{J}_{cdab})$
$\displaystyle+\frac{1}{12}\sum_{abcde}\sum_{J_{1}J_{2}\mathcal{J}}(n_{a}n_{b}\bar{n}_{c}\bar{n}_{d}\bar{n}_{e}+\bar{n}_{a}\bar{n}_{b}n_{c}n_{d}n_{e})\frac{\hat{\mathcal{J}}^{2}}{\hat{j}_{i}^{2}}(A^{J_{1}J_{2}\mathcal{J}}_{abicde}B^{J_{2}J_{1}\mathcal{J}}_{cdeabj}-B^{J_{1}J_{2}\mathcal{J}}_{abicde}A^{J_{2}J_{1}\mathcal{J}}_{cdeabj})$
$\displaystyle C^{J}_{ijkl}=$
$\displaystyle\sum_{a}(A_{ia}B^{J}_{ajkl}+A_{ja}B^{J}_{iakl}-B^{J}_{ijal}A_{ak}-B^{J}_{ijka}A_{al})$
(56)
$\displaystyle-\sum_{a}(B_{ia}A^{J}_{ajkl}+B_{ja}A^{J}_{iakl}-A^{J}_{ijal}B_{ak}-A^{J}_{ijka}B_{al})$
$\displaystyle+\frac{1}{2}\sum_{ab}\sum_{J}(\bar{n}_{a}\bar{n}_{b}-n_{a}n_{b})(A^{J}_{ijab}B^{J}_{abkl}-B^{J}_{ijab}A^{J}_{abkl})$
$\displaystyle+\sum_{ab}\sum_{J^{\prime}}(n_{a}-n_{b})\hat{J}^{{}^{\prime}2}(1-P^{J}_{ij})\begin{Bmatrix}i&j&J\\\
k&l&J^{\prime}\end{Bmatrix}(\bar{A}^{J^{\prime}}_{i\bar{l}a\bar{b}}\bar{B}^{J^{\prime}}_{a\bar{b}k\bar{j}}-\bar{B}^{J^{\prime}}_{i\bar{l}a\bar{b}}\bar{A}^{J^{\prime}}_{a\bar{b}k\bar{j}})$
$\displaystyle+\sum_{ab}\sum_{\mathcal{J}}(n_{a}-n_{b})\frac{\hat{\mathcal{J}}^{2}}{\hat{J}^{2}}(A_{ab}B^{JJ\mathcal{J}}_{ijbkla}-B_{ab}A^{JJ\mathcal{J}}_{ijbkla})$
$\displaystyle-\frac{1}{2}\sum_{abc}\sum_{J^{\prime}\mathcal{J}}(n_{a}n_{b}\bar{n}_{c}+\bar{n}_{a}\bar{n}_{b}n_{c})\frac{\hat{J}^{\prime}\hat{\mathcal{J}}^{2}}{\hat{J}}\Bigl{[}(1-P_{ij}^{J})\begin{Bmatrix}i&j&J\\\
c&\mathcal{J}&J^{\prime}\end{Bmatrix}(A^{J^{\prime}}_{cjab}B^{J^{\prime}J\mathcal{J}}_{abiklc}-B^{J^{\prime}}_{cjab}A^{J^{\prime}J\mathcal{J}}_{abiklc})$
$\displaystyle\hskip 150.00023pt-(1-P_{kl}^{J})\begin{Bmatrix}k&l&J\\\
c&\mathcal{J}&J^{\prime}\end{Bmatrix}(B^{JJ^{\prime}\mathcal{J}}_{ijcabk}A^{J^{\prime}}_{abcl}-A^{JJ^{\prime}\mathcal{J}}_{ijcabk}B^{J^{\prime}}_{abcl})\Bigr{]}$
$\displaystyle+\frac{1}{6}\sum_{abcd}\sum_{J^{\prime}\mathcal{J}}(n_{a}n_{b}n_{c}\bar{n}_{d}-\bar{n}_{a}\bar{n}_{b}\bar{n}_{c}n_{d})\frac{\hat{\mathcal{J}}^{2}}{\hat{J}^{2}}(A^{JJ^{\prime}\mathcal{J}}_{ijdabc}B^{J^{\prime}J\mathcal{J}}_{abckld}-B^{JJ^{\prime}\mathcal{J}}_{ijdabc}A^{J^{\prime}J\mathcal{J}}_{abckld})$
$\displaystyle+\frac{1}{4}\sum_{\begin{subarray}{c}abcd\\\
J_{ab}J_{cd}\end{subarray}}\sum_{J_{ph}}(\bar{n}_{a}\bar{n}_{b}n_{c}n_{d}-n_{a}n_{b}\bar{n}_{c}\bar{n}_{d})(1-P^{J}_{ij})(1-P^{J}_{kl})\hat{J_{ph}}^{2}\begin{Bmatrix}i&j&J\\\
k&l&J_{ph}\end{Bmatrix}\bar{A}^{J_{ph}}_{i\bar{l};(cdJ_{cd})(\overline{abJ_{ab}})}\bar{B}^{J_{ph}}_{(cdJ_{cd})(\overline{abJ_{ab}});k\bar{j}}$
$\displaystyle C^{J_{1}J_{2}\mathcal{J}}_{ijklmn}=$
$\displaystyle\sum_{a}\left\\{P^{J_{1}\mathcal{J}}_{ij/k}(A_{ka}B^{J_{1}J_{2}\mathcal{J}}_{ijalmn}-B_{ka}A^{J_{1}J_{2}\mathcal{J}}_{ijalmn})-P^{J_{2}\mathcal{J}}_{lm/n}(B^{J_{1}J_{2}\mathcal{J}}_{ijklma}A_{an}-A^{J_{1}J_{2}\mathcal{J}}_{ijklma}B_{an})\right\\}$
(57)
$\displaystyle+\sum_{a}P^{J_{1}\mathcal{J}}_{ij/k}P^{J_{2}\mathcal{J}}_{lm/n}\hat{J}_{1}\hat{J}_{2}\begin{Bmatrix}n&a&J_{1}\\\
k&\mathcal{J}&J_{2}\end{Bmatrix}(A^{J_{1}}_{ijna}B^{J_{2}}_{kalm}-B^{J_{1}}_{ijna}A^{J_{2}}_{kalm})$
$\displaystyle+\frac{1}{2}\sum_{ab}(\bar{n}_{a}\bar{n}_{b}-n_{a}n_{b})P^{J_{1}\mathcal{J}}_{ij/k}(A^{J_{1}}_{ijab}B^{J_{1}J_{2}\mathcal{J}}_{abklmn}-B^{J_{1}}_{ijab}A^{J_{1}J_{2}\mathcal{J}}_{abklmn})$
$\displaystyle-\frac{1}{2}\sum_{ab}(\bar{n}_{a}\bar{n}_{b}-n_{a}n_{b})P^{J_{2}\mathcal{J}}_{lm/n}(B^{J_{1}J_{2}\mathcal{J}}_{ijkabn}A^{J_{2}}_{ablm}-A^{J_{1}J_{2}\mathcal{J}}_{ijkabn}B^{J_{2}}_{ablm})$
$\displaystyle-\sum_{ab}\sum_{J^{\prime}\mathcal{J}^{\prime}}(n_{a}\bar{n}_{b}-\bar{n}_{a}n_{b})\hat{J}^{{}^{\prime}2}\hat{\mathcal{J}}^{{}^{\prime}2}P^{J_{1}\mathcal{J}}_{ij/k}P^{J_{2}\mathcal{J}}_{lm/n}(-1)^{k+n+J_{1}+J_{2}+\mathcal{J}}\begin{Bmatrix}b&\mathcal{J}^{\prime}&J_{2}\\\
k&J_{1}&\mathcal{J}\\\
J^{\prime}&a&n\end{Bmatrix}(A^{J^{\prime}}_{bkan}B^{J_{1}J_{2}\mathcal{J}^{\prime}}_{ijalmb}-B^{J^{\prime}}_{bkan}A^{J_{1}J_{2}\mathcal{J}^{\prime}}_{ijalmb})$
$\displaystyle+\frac{1}{6}\sum_{abc}\sum_{J^{\prime}}(n_{a}n_{b}n_{c}+\bar{n}_{a}\bar{n}_{b}\bar{n}_{c})(A^{J_{1}J^{\prime}\mathcal{J}}_{ijkabc}B^{J^{\prime}J_{2}\mathcal{J}}_{abclmn}-B^{J_{1}J^{\prime}\mathcal{J}}_{ijkabc}A^{J^{\prime}J_{2}\mathcal{J}}_{abclmn})$
$\displaystyle+\frac{1}{2}\sum_{abc}\sum_{J^{\prime}\mathcal{J}^{\prime}\mathcal{J}^{\prime\prime}}(n_{a}n_{b}\bar{n}_{c}+\bar{n}_{a}\bar{n}_{b}n_{c})\hat{\mathcal{J}}^{{}^{\prime}2}\hat{\mathcal{J}}^{{}^{\prime\prime}2}P^{J_{1}\mathcal{J}}_{ij/k}P^{J_{2}\mathcal{J}}_{lm/n}\begin{Bmatrix}k&J^{\prime}&\mathcal{J}^{\prime}\\\
J_{1}&\mathcal{J}^{\prime\prime}&c\\\
\mathcal{J}&n&J_{2}\end{Bmatrix}(A^{J^{\prime}J_{2}\mathcal{J}^{\prime}}_{abklmc}B^{J_{1}J^{\prime}\mathcal{J}^{\prime\prime}}_{ijcabn}-B^{J^{\prime}J_{2}\mathcal{J}^{\prime}}_{abklmc}A^{J_{1}J^{\prime}\mathcal{J}^{\prime\prime}}_{ijcabn})$

##  Appendix C Normal-ordered commutator for the LGM model

Operators built from the usual quasi-spin operators $K_{\pm},K_{z}$ can be re-
expressed in terms of operators normal-ordered with respect to the mean-field
ground state in the unbroken phase (we indicate the normal-ordered operators
with braces $\\{\\}$) :

$\displaystyle K_{z}$ $\displaystyle=\\{K_{z}\\}-\tfrac{1}{2}N$ (58)
$\displaystyle K_{z}^{2}$
$\displaystyle=\\{K_{z}^{2}\\}+(\tfrac{1}{2}\\!-\\!N)\\{K_{z}\\}+\tfrac{1}{4}N^{2}$
(59) $\displaystyle K_{z}^{3}$
$\displaystyle=\\{K_{z}^{3}\\}+\tfrac{3}{2}(1-N)\\{K_{z}^{2}\\}$
$\displaystyle~{}~{}+\tfrac{1}{4}(1-3N+3N^{2})\\{K_{z}\\}-\tfrac{1}{8}N^{3}$
(60) $\displaystyle K_{\pm}^{n}$ $\displaystyle=\\{K_{\pm}^{n}\\}$ (61)
$\displaystyle K_{+}^{2}K_{z}$
$\displaystyle=\\{K_{+}^{2}K_{z}\\}-\tfrac{1}{2}N\\{K_{+}^{2}\\}$ (62)
$\displaystyle K_{z}K_{-}^{2}$
$\displaystyle=\\{K_{z}K_{-}^{2}\\}-\tfrac{1}{2}N\\{K_{-}^{2}\\}$ (63)
$\displaystyle K_{+}K_{-}$ $\displaystyle=\\{K_{+}K_{-}\\}$ (64)
$\displaystyle K_{-}K_{+}$ $\displaystyle=\\{K_{+}K_{-}\\}-2\\{K_{z}\\}+N$
(65) $\displaystyle K_{+}K_{z}K_{-}$
$\displaystyle=\\{K_{+}K_{z}K_{-}\\}-\tfrac{1}{2}N\\{K_{+}K_{-}\\}$ (66)

Relevant commutator expression in terms of the vacuum normal ordered operators
are

$\displaystyle[K_{z},K_{\pm}]$ $\displaystyle=\pm K_{\pm}$ (67)
$\displaystyle[K_{+},K_{-}]$ $\displaystyle=2K_{z}$ (68)
$\displaystyle[K_{+}^{2},K_{-}^{2}]$
$\displaystyle=8K_{+}K_{z}K_{-}+8K_{+}K_{-}$
$\displaystyle~{}~{}-8K_{z}^{2}-4K_{z}$ (69)
$\displaystyle[K_{z},K_{\pm}^{2}]$ $\displaystyle=\pm 2K_{\pm}^{2}$ (70)
$\displaystyle[K_{z}^{2},K_{+}^{2}]$
$\displaystyle=4K_{+}^{2}K_{z}+4K_{+}^{2}$ (71)
$\displaystyle[K_{z}^{2},K_{-}^{2}]$
$\displaystyle=-4K_{z}K_{-}^{2}-4K_{-}^{2}$ (72)
$\displaystyle[K_{+}K_{-},K_{+}^{2}]$
$\displaystyle=-4K_{+}^{2}K_{z}-2K_{+}^{2}$ (73)
$\displaystyle[K_{+}K_{-},K_{-}^{2}]$
$\displaystyle=4K_{z}K_{-}^{2}+2K_{-}^{2}$ (74)

In terms of the operators normal ordered with respect to the mean field ground
state, the commutator expressions are

$\displaystyle[\\{K_{z}\\},\\{K_{\pm}^{2}\\}]=$ $\displaystyle\pm
2\\{K_{\pm}^{2}\\}$ (76) $\displaystyle[\\{K_{z}^{2}\\},\\{K_{\pm}^{2}\\}]=$
$\displaystyle\pm 4\\{K_{\pm}^{2}K_{z}\\}\pm 3\\{K_{\pm}^{2}\\}$ (77)
$\displaystyle[\\{K_{+}K_{-}\\},\\{K_{\pm}^{2}\\}]=$ $\displaystyle\mp
4\\{K_{\pm}^{2}K_{z}\\}$ $\displaystyle~{}\pm 2(N\\!-\\!1)\\{K_{\pm}^{2}\\}$
(78) $\displaystyle[\\{K_{+}^{2}\\},\\{K_{-}^{2}\\}]=$ $\displaystyle
8\\{K_{+}K_{z}K_{-}\\}$
$\displaystyle~{}+(8\\!-\\!4N)\\{K_{+}K_{-}\\}-8\\{K_{z}^{2}\\}$
$\displaystyle~{}+8(N\\!-\\!1)\\{K_{z}\\}$ $\displaystyle~{}-2N(N\\!-\\!1)$
(79) $\displaystyle[\\{K_{+}K_{z}K_{-}\\},\\{K_{+}^{2}\\}]$
$\displaystyle=2\\{K_{+}^{3}K_{-}\\}-4\\{K_{+}^{2}K_{z}^{2}\\}$
$\displaystyle~{}+2(N\\!-\\!4)\\{K_{+}^{2}K_{z}\\}$
$\displaystyle~{}+2(N\\!-\\!1)\\{K_{+}^{2}\\}.$ (80)

##  Appendix D CCD solution to the LGM model

The coupled cluster doubles solution of the LGM model is obtained from the
decoupling equation

$\langle\phi|\\{K_{-}^{2}\\}e^{-T}He^{T}|\phi\rangle=0.$ (81)

where the CCD excitation operator is $T=t_{2}\\{K_{+}^{2}\\}$. To solve this,
we need the coefficient of $\\{K_{+}^{2}\\}$ in $H$, $[H,T]$, and $[[H,T],T]$.
Using the expressions in (67), we obtain

$\tfrac{1}{2}V+2\epsilon t_{2}+2V(N^{2}-7N+9)t_{2}^{2}=0$ (82)

defining $\alpha\equiv N^{2}-7N+9$, we obtain the solution

$t_{2}=\frac{-\epsilon}{2V\alpha}\left[1-\sqrt{1-\alpha\left(\frac{V}{\epsilon}\right)^{2}}\right].$
(83)

We are interested in the region where $NV/\epsilon\sim 1$, so
$t_{2}\sim\frac{N}{\alpha}\to\frac{1}{N}$ as $N\to\infty$. The CCD correlation
energy is

$\displaystyle E_{CCD}$ $\displaystyle=\langle\phi|[H,T]\phi\rangle$ (84)
$\displaystyle=\tfrac{1}{2}Vt_{2}\langle\phi|[\\{K_{-}^{2}\\},\\{K_{+}^{2}\\}]|\phi\rangle$
$\displaystyle=Vt_{2}N(N-1)$

which goes like $N$ as $N\to\infty$.

##  Appendix E Magnus(2) solution to the LGM model

In the Magnus formulation of the IMSRG, the transformed Hamiltonian is given
by (8). At the NO2B level, the only off-diagonal operators that respect the
symmetries of the problem are $K_{\pm}^{2}$. The requirement of
antihermiticity means $\eta(s)\propto(K_{+}^{2}-K_{-}^{2})$ for all $s$, and
so all the commutators in (6) vanish. Consequently, $\Omega$ is also
proportional to this operator, and the transformation is specified by a single
real number:

$\Omega=\tfrac{1}{2}\Omega_{\rm 2b}(K_{+}^{2}-K_{-}^{2}).$ (85)

Next we consider the behavior of nested commutators of $\Omega$ with $H$.
$[\Omega,H]$ yields one term proportional to $K_{+}^{2}+K_{-}^{2}$, and a
second term proportional to

$[K_{+}^{2}-K_{-}^{2},K_{+}^{2}+K_{-}^{2}]=2[K_{+}^{2},K_{-}^{2}],$ (86)

so that

$[\Omega,H]=a_{1}[K_{+}^{2},K_{-}^{2}]+b_{1}(K_{+}^{2}+K_{-}^{2}).$ (87)

By applying the commutator expressions in (67) we find in the NO2B
approximation the double commutator

$[K_{\pm}^{2},[K_{+}^{2},K_{-}^{2}]]=8\tilde{\alpha}K_{\pm}^{2}$ (88)

with $\tilde{\alpha}\equiv N^{2}-5N+7$. From (87), (86) and (88), we see that
arbitrary nested commutators $[\Omega,\ldots[\Omega,H]\ldots]$ will have the
same structure as (87), so we may write for the $k$-fold nested commutator

$[\Omega,H]^{(k)}=a_{k}[K_{+}^{2},K_{-}^{2}]+b_{k}(K_{+}^{2}+K_{-}^{2}).$ (89)

We have a recurrence relation for $a_{k}$ and $b_{k}$

$a_{k}=\theta^{2}a_{k-2}~{},\hskip 20.00003ptb_{k}=\theta^{2}b_{k-2}~{}$ (90)

with $\theta\equiv\sqrt{4\Omega_{\rm 2b}^{2}\tilde{\alpha}}$, and initial
values

$\displaystyle a_{1}$ $\displaystyle=\tfrac{1}{2}\Omega_{\rm 2b}V\hskip
10.00002pt$ $\displaystyle a_{2}=-\Omega_{\rm 2b}^{2}\epsilon$ (91)
$\displaystyle b_{0}$ $\displaystyle=\tfrac{1}{2}V$ $\displaystyle
b_{1}=-\Omega_{\rm 2b}\epsilon.$ (92)

Writing the off-diagonal piece of the transformed Hamiltonian as
$\tilde{H}^{\rm od}=\tfrac{1}{2}\tilde{V}(K_{+}^{2}+K_{-}^{2})$, we have

$\begin{aligned}
\tfrac{1}{2}\tilde{V}&=\sum_{k=0}^{\infty}b_{k}\frac{1}{k!}\\\
&=\tfrac{1}{2}V\cosh\theta-\Omega_{\rm
2b}\epsilon\frac{\sinh\theta}{\theta}\end{aligned}.$ (93)

The requirement that $\tilde{H}^{\rm od}\to 0$ gives the decoupling condition

$\tanh\theta=\sqrt{\tilde{\alpha}}\frac{V}{\epsilon}$ (94)

which fixes $\Omega_{\rm 2b}$. The ground state energy is obtained from the
zero-body piece of $\tilde{H}$. From the commutator expressions (67) the zero-
body piece of all nested commutator comes from $[K_{+}^{2},K_{-}^{2}]_{\rm
0b}=-2N(N-1)$. Consequently, the IMSRG(2) correlation energy is

$\displaystyle E_{\rm IMSRG(2)}$
$\displaystyle=-2N(N-1)\left(\frac{V}{4\sqrt{\tilde{\alpha}}}\sinh\theta+\frac{\epsilon}{4\tilde{\alpha}}(1-\cosh\theta)\right)$
(95)
$\displaystyle=-\frac{N(N-1)}{2\tilde{\alpha}}\left(1-\sqrt{1-\tilde{\alpha}\frac{V^{2}}{\epsilon^{2}}}\right).$

In the second line of (95) we have used (94) and some identities of hyperbolic
functions.

## References

  * Hergert [2020] H. Hergert, A guided tour of ab initio nuclear many-body theory, Frontiers in Physics 8, 379 (2020).
  * Miyagi  _et al._ [2022] T. Miyagi, S. R. Stroberg, P. Navrátil, K. Hebeler, and J. D. Holt, Converged ab initio calculations of heavy nuclei, Physical Review C 105, 014302 (2022).
  * Hu  _et al._ [2022] B. Hu, W. Jiang, T. Miyagi, Z. Sun, A. Ekström, C. Forssén, G. Hagen, J. D. Holt, T. Papenbrock, S. R. Stroberg, _et al._ , Ab initio predictions link the neutron skin of 208pb to nuclear forces, Nature Physics 18, 1196 (2022).
  * Hebeler  _et al._ [2023] K. Hebeler, V. Durant, J. Hoppe, M. Heinz, A. Schwenk, J. Simonis, and A. Tichai, Normal ordering of three-nucleon interactions for ab initio calculations of heavy nuclei, Physical Review C 107, 024310 (2023).
  * Hergert  _et al._ [2013] H. Hergert, S. K. Bogner, S. Binder, A. Calci, J. Langhammer, R. Roth, and A. Schwenk, In-medium similarity renormalization group with chiral two- plus three-nucleon interactions, Phys. Rev. C 87, 034307 (2013).
  * Somà  _et al._ [2014] V. Somà, A. Cipollone, C. Barbieri, P. Navrátil, and T. Duguet, Chiral two- and three-nucleon forces along medium-mass isotope chains, Physical Review C 89, 061301(R) (2014).
  * Sun  _et al._ [2018] Z. H. Sun, T. D. Morris, G. Hagen, G. R. Jansen, and T. Papenbrock, Shell-model coupled-cluster method for open-shell nuclei, Physical Review C 98, 054320 (2018).
  * Stroberg  _et al._ [2021] S. R. Stroberg, J. D. Holt, A. Schwenk, and J. Simonis, Ab initio limits of atomic nuclei, Physical Review Letters 126, 022501 (2021).
  * Tichai  _et al._ [2023] A. Tichai, P. Demol, and T. Duguet, Towards heavy-mass ab initio nuclear structure: Open-shell ca, ni and sn isotopes from bogoliubov coupled-cluster theory, arXiv preprint arXiv:2307.15619 https://doi.org/10.48550/arXiv.2307.15619 (2023).
  * Frame  _et al._ [2018] D. Frame, R. He, I. Ipsen, D. Lee, D. Lee, and E. Rrapaj, Eigenvector continuation with subspace learning, Physical review letters 121, 032501 (2018).
  * Ekström and Hagen [2019] A. Ekström and G. Hagen, Global sensitivity analysis of bulk properties of an atomic nucleus, Physical Review Letters 123, 252501 (2019).
  * Melendez  _et al._ [2019] J. A. Melendez, R. J. Furnstahl, D. R. Phillips, M. T. Pratola, and S. Wesolowski, Quantifying correlated truncation errors in effective field theory, Physical Review C 100, 044001 (2019).
  * Jiang and Forssén [2022] W. Jiang and C. Forssén, Bayesian probability updates using sampling/importance resampling: Applications in nuclear theory, Frontiers in Physics 10, 1058809 (2022).
  * Belley  _et al._ [2024] A. Belley, J. M. Yao, B. Bally, J. Pitcher, J. Engel, H. Hergert, J. D. Holt, T. Miyagi, T. R. Rodríguez, A. M. Romero, S. R. Stroberg, and X. Zhang, Ab initio uncertainty quantification of neutrinoless double-beta decay in ${}^{76}\mathrm{Ge}$, Phys. Rev. Lett. 132, 182502 (2024).
  * Tsukiyama  _et al._ [2011] K. Tsukiyama, S. K. Bogner, and A. Schwenk, In-Medium Similarity Renormalization Group For Nuclei, Phys. Rev. Lett. 106, 222502 (2011).
  * Tsukiyama  _et al._ [2012] K. Tsukiyama, S. K. Bogner, and A. Schwenk, In-medium similarity renormalization group for open-shell nuclei, Phys. Rev. C 85, 061304(R) (2012).
  * Hergert  _et al._ [2016] H. Hergert, S. K. Bogner, T. D. Morris, A. Schwenk, and K. Tsukiyama, The In-Medium Similarity Renormalization Group: A Novel Ab Initio Method for Nuclei, Phys. Rep. 621, 165 (2016).
  * Hergert [2016] H. Hergert, In-medium similarity renormalization group for closed and open-shell nuclei, Physica Scripta 92, 023002 (2016).
  * Stroberg  _et al._ [2019] S. R. Stroberg, S. K. Bogner, H. Hergert, and J. D. Holt, Nonempirical Interactions for the Nuclear Shell Model: An Update, Annu. Rev. Nucl. Part. Sci. 69, 307 (2019), arXiv:1902.06154 .
  * White [2002] S. R. White, Numerical canonical transformation approach to quantum many-body problems, J. Chem. Phys. 117, 7472 (2002).
  * Hergert  _et al._ [2017] H. Hergert, S. K. Bogner, J. G. Lietz, T. D. Morris, S. J. Novario, N. M. Parzuchowski, and F. Yuan, In-Medium Similarity Renormalization Group Approach to the Nuclear Many-Body Problem, in  _Lect. Notes Phys._, edited by M. Hjorth-Jensen, M. P. Lombardo, and U. van Kolck (Springer, 2017) pp. 477–570, arXiv:1612.08315 .
  * Morris  _et al._ [2015] T. D. Morris, N. M. Parzuchowski, and S. K. Bogner, Magnus expansion and in-medium similarity renormalization group, Phys. Rev. C 92, 034331 (2015).
  * Jurgenson  _et al._ [2009] E. D. Jurgenson, P. Navrátil, and R. J. Furnstahl, Evolution of Nuclear Many-Body Forces with the Similarity Renormalization Group, Phys. Rev. Lett. 103, 082501 (2009).
  * Roth  _et al._ [2014] R. Roth, A. Calci, J. Langhammer, and S. Binder, Evolved chiral NN+3N Hamiltonians for ab initio nuclear structure calculations, Phys. Rev. C 90, 024325 (2014).
  * Hagen  _et al._ [2014] G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, Coupled-cluster computations of atomic nuclei, Reports on Progress in Physics 77, 096302 (2014).
  * Roth  _et al._ [2012] R. Roth, S. Binder, K. Vobig, A. Calci, J. Langhammer, and P. Navrátil, Medium-mass nuclei with normal-ordered chiral NN+3N interactions, Physical Review Letters 109, 052501 (2012).
  * Djärv  _et al._ [2021] T. Djärv, A. Ekström, C. Forssén, and G. R. Jansen, Normal-ordering approximations and translational (non) invariance, Physical Review C 104, 024324 (2021).
  * Kehrein [2006] S. Kehrein,  _The flow equation approach to many-particle systems_ , Vol. 217 (Springer Science & Business Media, 2006).
  * Machleidt and Entem [2011] R. Machleidt and D. R. Entem, Chiral effective field theory and nuclear forces, Physics Reports 503, 1 (2011).
  * Epelbaum  _et al._ [2009] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Modern theory of nuclear forces, Reviews of Modern Physics 81, 1773 (2009).
  * Kriel  _et al._ [2005] J. N. Kriel, A. Y. Morozov, and F. G. Scholtz, Non-perturbative flow equations from continuous unitary transformations, J. Phys. A. Math. Gen. 38, 205 (2005), arXiv:0602022 [cond-mat] .
  * Friman and Schwenk [2011] B. Friman and A. Schwenk, Three-body interactions in fermi systems, in  _From Nuclei To Stars: Festschrift in Honor of Gerald E Brown_ (World Scientific, 2011) pp. 141–156.
  * Tsukiyama [2010] K. Tsukiyama,  _Similarity Renormalization Group for Nuclear Many-Body Systems_ , Ph.D. thesis, The University of Tokyo (2010).
  * Morris [2016] T. D. Morris,  _A Novel Approach to Solving the In-Medium SRG_ , Ph.D. thesis, Michigan State University (2016).
  * Heinz  _et al._ [2021] M. Heinz, A. Tichai, J. Hoppe, K. Hebeler, and A. Schwenk, In-medium similarity renormalization group with three-body operators, Phys. Rev. C 103, 044318 (2021).
  * Shavitt and Bartlett [2009] I. Shavitt and R. J. Bartlett,  _Many-Body Methods in Chemistry and Physics_ (Cambridge University Press, New York, 2009).
  * Evangelista and Gauss [2012] F. A. Evangelista and J. Gauss, On the approximation of the similarity-transformed hamiltonian in single-reference and multireference coupled cluster theory, Chemical Physics 401, 27 (2012).
  * Evangelista [2014] F. A. Evangelista, A driven similarity renormalization group approach to quantum many-body problems, The Journal of Chemical Physics 141, 054109 (2014).
  * Li and Evangelista [2020] C. Li and F. A. Evangelista, Connected three-body terms in single-reference unitary many-body theories: Iterative and perturbative approximations, The Journal of Chemical Physics 152, 234116 (2020).
  * He and Stroberg [2024] B. C. He and S. R. Stroberg, Factorized approximation to the imsrg(3), arXiv preprint arXiv:2405.19594 https://doi.org/10.48550/arXiv.2405.19594 (2024).
  * [41] S. R. Stroberg, IMSRG with flowing 3N operators, https://abinitio.triumf.ca/2020/stroberg_AbInitio2020.pdf.
  * Lipkin  _et al._ [1965] H. J. Lipkin, N. Meshkov, and A. J. Glick, Validity of many-body approximation methods for a solvable model. (i). exact solutions and perturbation theory, Nucl. Phys. 62, 188 (1965).
  * Lührmann [1977] K. Lührmann, Equations for subsystems, Annals of Physics 103, 253 (1977).
  * Harsha  _et al._ [2018] G. Harsha, T. Shiozaki, and G. E. Scuseria, On the difference between variational and unitary coupled cluster theories, The Journal of Chemical Physics 148, 10.1063/1.5011033 (2018).
  * Pirner and Friman [1998] H. J. Pirner and B. Friman, Hamiltonian flow equations for the Lipkin model, Phys. Lett. Sect. B Nucl. Elem. Part. High-Energy Phys. 434, 231 (1998).
  * Mielke [1998] A. Mielke, Flow equations for band-matrices, The European Physical Journal B-Condensed Matter and Complex Systems 5, 605 (1998).
  * Kriel  _et al._ [2007] J. N. Kriel, F. G. Scholtz, and J. D. Thom, Moyal implementation of flow equations - A non-perturbative approach to quantum many-body systems, J. Phys. A Math. Theor. 40, 9483 (2007).
  * Scholtz  _et al._ [2003] F. G. Scholtz, B. H. Bartlett, and H. B. Geyer, Nonperturbative flow equations from running expectation values, Phys. Rev. Lett. 91, 080602 (2003).
  * Dusuel and Vidal [2005] S. Dusuel and J. Vidal, Continuous unitary transformations and finite-size scaling exponents in the Lipkin-Meshkov-Glick model, Phys. Rev. B 71, 224420 (2005).
  * Neuscamman  _et al._ [2009] E. Neuscamman, T. Yanai, and G. K.-L. Chan, Quadratic canonical transformation theory and higher order density matrices, The Journal of Chemical Physics 130, 124102 (2009).
  * Hagen  _et al._ [2009] G. Hagen, T. Papenbrock, and D. J. Dean, Solution of the center-of-mass problem in nuclear structure calculations, Physical Review Letters 103, 062503 (2009).
  * Tichai  _et al._ [2020] A. Tichai, R. Wirth, J. Ripoche, and T. Duguet, Symmetry reduction of tensor networks in many-body theory: I. automated symbolic evaluation of su (2) algebra, The European Physical Journal A 56, 1 (2020).
  * Stroberg [2018] S. R. Stroberg, IMSRG++ code (2018).
  * Sanderson and Curtin [2016] C. Sanderson and R. Curtin, Armadillo: a template-based C++ library for linear algebra, Journal of Open Source Software 1, 26 (2016).
  * Sanderson and Curtin [2019] C. Sanderson and R. Curtin, Practical Sparse Matrices in C++ with Hybrid Storage and Template-Based Expression Optimisation, Mathematical and Computational Applications 24, 10.3390/mca24030070 (2019).