EFFECTIVE BETA–DECAY OPERATORS WITH COUPLED CLUSTER
THEORY
By
Samuel John Novario

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics—Doctor of Philosophy
2018

ABSTRACT
EFFECTIVE BETA–DECAY OPERATORS WITH COUPLED CLUSTER THEORY
By
Samuel John Novario
Coupled Cluster theory is a powerful ab initio framework for solving the many-body the
Schrödinger equation and has been utilized successfully to describe the highly-correlated
systems found in quantum chemistry and nuclear physics. This method uses a special similarity transformation to decouple a system’s ground state from excitations from it. This
transformation contains significant correlations that can be used to extend coupled cluster
theory to excited states and open-shell systems with the equation-of-motions method. Additionally, properties of these states can be obtained by consistently transforming relevant
operators using the coupled cluster similarity transformation. The coupled cluster method is
systematically improvable and scales polynomially with the system size. With this flexibility
and reach, coupled cluster theory can be applied across the nuclear chart to contribute to
many important open problems in physics.
Several fundamental questions in modern physics involve electroweak interactions within
nuclei, including the search for the elusive neutrinoless double-beta decay. Often the largest
uncertainty within these experiments is due to nuclear-structure-dependent quantities that
are calculated within some many-body framework. The main focus of this thesis is to apply the coupled cluster method to calculate effective Fermi and Gamow-Teller beta-decay
operators between open shell states. By confirming the validity of this method, it can be
extended to double-beta decay and other electroweak processes.

ACKNOWLEDGMENTS

First and foremost, I would like to thank my parents, Amy and Steve, for their unconditional support of my educational aspirations and passionate cultivation of my curiosities in
addition to their persistent emotional, financial, and logistical aid they have provided during
my long journey up till this point. Also, I would like to thank my siblings, Maggie, Lucy,
and Charlie, for their friendship and inspiration, as well as the examples they set which I was
able to follow. Moreover, I would like to thank my grandparents, aunts, uncles and cousins
for their words of encouragement and praise. Lastly, I would like to thank my friends in East
Lansing, Illinois, and around the country for helping me have fun and maintain my peace of
mind, and never failing to ask me when I would graduate.
Next, I would like to thank my all my teachers throughout my long tenure as a student.
I would especially like to acknowledge my science and math teachers that fostered and
encouraged my interests. In particular, I would like to credit Bryon Leonard for opening
my mind to the beauty of modern physics and the path that leads to where I am now.
Additionally, I am thankful for Philippe Collon and Michael Thoennessen for giving me my
first research opportunities as an undergraduate student when I’m sure I was more of a
burden than a blessing.
I am grateful to Dr. Georg Bollen and Dr. Ryan Ringle for inviting me into the LEBIT
collaboration, my first academic family. Along with my fellow graduate students, David
Lincoln, Scott Bustabad, and Brad Barquest, my LEBIT colleagues helped me find my place
as a graduate student, guided me as a scientist and collaborator, provided opportunities to
travel and present at national conferences, and gave me research experience that few theorists
hold. While in LEBIT, I was fortunate to share an office with Jenna Smith, who looked over

iii

me as a young graduate student and gave me valuable advice and encouragement when it
was greatly needed.
When I mentioned that I decided to change my research direction, my current advisor,
Morten Hjorth-Jensen, without hesitation and without any solicitation from me, took a risk
by taking me on as a student. For this, I am eternally grateful, and because Morten has
become an invaluable mentor, colleague, and friend, I feel very lucky for finding myself in
that situation. Along with my advisor, I would like to thank the members of my committee–
Scott Bogner, Piotr Piecuch, Carlo Piermarocchi, and Remco Zegers–for agreeing to join my
adventure through graduate school and guiding me along the way. I’m also indebted to Heiko
Hergert for fielding my relentless questions, Ragnar Stroberg for unknowingly helping with
his clean and well-commented code, and Kim Crosslan for shouldering most of my blunders
for more than 7 years and for making sure I was aware of them.
Lasly, I would like to thank my fellow graduate students. To Justin Lietz, my office mate,
world traveling partner, coupled-cluster colleague, and friend, I’m thankful for lending his
advice, his expertise, and his ear. To Nathan Parzchuowki, Titus Morris, and Fei Yuan, my
fellow nuclear theory colleagues and friends, I’m thankful for their constant help and the
pleasant work environment that they created.
And to all of the people that I haven’t mentioned, I realize that I am just a product of
the circumstances that unfolded around me, especially the people with which I’ve crossed
paths. So I’d like to give thanks to everyone, and to no one in particular, who helped shape
my success in completing this dissertation.

iv

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . . . . . . . . . . . . xiii
Chapter 1 Introduction . . . . . . . . . . . . . .
1.1 A Brief History of Nuclear Structure Theory
1.2 Electroweak Theory and Nuclear Structure .
1.3 Ab-Initio Descriptions of Beta Decay . . . .
1.4 Thesis Structure . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

1
2
4
6
7

Chapter 2 Many-Body Quantum Mechanics
2.1 Independent-Particle Model . . . . . . . .
2.2 Second Quantization . . . . . . . . . . . .
2.3 Normal Ordering . . . . . . . . . . . . . .
2.4 Wick’s Theorem . . . . . . . . . . . . . . .
2.5 Hartree–Fock Method . . . . . . . . . . . .
2.6 Configuration-Interaction . . . . . . . . . .
2.7 Many-Body Perturbation Theory . . . . .
2.7.1 Factorization Theorem . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

8
8
10
12
15
20
25
27
31

Chapter 3 Coupled-Cluster Theory . . . . . .
3.1 Exponential Ansatz . . . . . . . . . . . . . .
3.1.1 The Coupled Cluster Equations . . .
3.2 Linked-Cluster Theorem and MBPT . . . .
3.3 Example: Pairing Model . . . . . . . . . . .
3.4 Solving the Coupled Cluster Equations . . .
3.4.1 Symmetry Channels . . . . . . . . .
3.4.2 Matrix Structures and Intermediates
3.5 Example: Homogeneous Electron Gas . . . .
3.6 Coupled Cluster for Finite Nuclei . . . . . .
3.6.1 Harmonic Oscillator Basis . . . . . .
3.6.2 The Nuclear Interaction . . . . . . .
3.6.3 Ground-State Results for Nuclei . . .
3.7 Ground-State Center-of-Mass Factorization .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

34
34
37
38
43
48
49
51
57
62
62
65
68
71

Chapter 4 Equation-of-Motion Method .
4.1 Equation-of-Motion States . . . . . . .
4.2 Dual Solutions . . . . . . . . . . . . . .
4.2.1 Induced Three-Body Interaction
4.3 Solving the EOM equations . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

76
77
80
82
83

v

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

4.4

4.5
4.6

EOM-CC for 2D Quantum Dots
4.4.1 Quantum Dot Formalism
4.4.2 Quantum Dot Results .
Quality of EOM Solutions . . .
EOM-CC for Finite Nuclei . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

90
90
92
93
96

Chapter 5 Beta–Decay Effective Operators . . . . . . . . . . . . . . . . . . . 103
5.1 Beta-Decay Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.2 Coupled Cluster Effective Operators . . . . . . . . . . . . . . . . . . . . . . 108
Chapter 6 Conclusions and Perspectives . . . . . . . . . . . . . . . . . . . . . 114
6.1 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
APPENDICES . . . . . . . . . . . . . . . . .
Appendix A CCSD Diagrams . . . . . . . .
Appendix B Effective Hamiltonian Diagrams
Appendix C Computational Implementation
Appendix D Angular Momentum Coupling .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

116
117
124
131
142

REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

vi

LIST OF TABLES

Table 5.1:

Summary of the selection rules for allowed beta decays according to the
angular momentum J and parity π of the initial (I) and final (F ) states.

vii

106

LIST OF FIGURES

Figure 1.1:

Figure 1.2:

Figure 2.1:

Figure 2.2:

Figure 2.3:

Figure 3.1:

Figure 3.2:

Figure 3.3:

Nuclear chart of nuclei with ground-state energies which have been calculated with ab-initio methods and NN+3N interactions. Figure taken
from [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

Progress of ab-initio nuclear structure from calculations of ground-state
energies with NN+3N interactions. Early progress was approximately
linear as the problem size scaled with Moore’s law while more recent
progress has taken advantage of new algorithms which have outpaced
Moore’s law. Data taken from [1]. . . . . . . . . . . . . . . . . . . . . .

5

A depiction of the closed-shell reference state in the independent particle
model. Each horizontal line represents a shell of single-particle orbits,
represented by circles, and the dotted line represents the Fermi level
which separates the unoccupied particle states from the occupied hole
states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

A depiction of 1p -1h, 2p -2h, 1p -0h, and 0p -1h Slater determinants defined relative to the reference state in the independent particle model. .

14

Scaling of the matrix size and number of non-zero matrix elements for
nuclear CI calculations of light nuclei. Even for modestly-sized model
spaces, the memory requirements approach the limit of petascale supercomputers (∟ 1010 ). Figure taken from [2]. . . . . . . . . . . . . . . . .

27

Schematic representation of the pairing model space. The shells are
equally spaced and doubly degenerate with one spin-up and one spindown state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

Correlation energy for the pairing model with exact diagonalization, CCD,
and perturbation theory to third (MBPT3) and fourth order (MBPT4)
for a range of interaction values, g. . . . . . . . . . . . . . . . . . . . . .

46

Visualization of the CCD similarity transform on the pairing Hamiltonian
for four particles and six shells. This shows the main function of CCD,
which is to decouple 2p -2h excitations from the ground state, shown by
the suppression of matrix elements on the first column. In the pairing
model, this also has the effect of decoupling 2p -2h excitations from 4p -4h
excitations. Also, the non-unitary nature of the transformation is obvious
given the asymmetry of the resulting Hamiltonian. . . . . . . . . . . . .

48

viii

Figure 3.4:

Figure 3.5:

Visulization of the Fourier transform of a finite box. This transformation
characterizes the construction of the single-particle basis for infinite matter, mapping plane waves in coordinate space onto finitely-spaced points
in momentum space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

CCD energy per electron in Hartrees for the 3D homogeneous electron
gas as function of the Wigner-Seitz radius in units of Bohr radii. The
calculation used periodic boundary conditions and a basis with 25 shells,
resulting in a total of 1238 single-particle states. Also plotted are the
variational quantum Monte Carlo (VMC) results from [3]. . . . . . . . .

61

Figure 3.6:

A schematic illustration of the harmonic oscillator basis used for calculations of nuclei. Shown is an example of a initial reference state for carbon14, with 6 protons filled to the p3/2-subshell closure and 8 neutrons filled
to the p1/2-shell closure. See text for details on the single-particle states. 64

Figure 3.7:

Diagrammatic form of the chiral EFT expansion up to N3 LO. The solid
lines represent nucleons and the dashed lines represent pions. The different vertices represent higher-order interactions. Figure taken from [4]. .

67

Ground-state energies for 16 O for the EM N3 LO NN only interaction
and with the added 3N interaction from Navrátil, both SRG softened
with λSRG = 1.88, 2.24 fm−1 . The energies are plotted for emax = 10, 12.
The most obvious difference is between the NN and NN+3N calculations,
showing the importance of including 3N forces. The differences between
the cutoff parameters are resolved within ∟ 1% with the inclusion of 3N
forces and can be rectified further by including additional correlations or
full 3N forces. The experimental binding energy is shown with the grey
dashed line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

Ground-state energies for doubly magic nuclei as a function of the harmonic oscillator energy ~ω with the NN+3N(400) interaction, SRG softened with λSRG = 2fm−1 . The energies are plotted for emax = 8, 10, 12,
showing the convergence as the model space increases. The results are independent of the underlying oscillator frequency to ∟ 1% for emax = 12.
The grey dashed line is the experimental binding energy. The overbinding
of this interaction becomes apparent as the system size increases. . . . .

70

Figure 3.10: Ground-state energies for singly magic nuclei as a function of the harmonic oscillator energy ~ω with the NN+3N(400) interaction, SRG softened with λSRG = 2fm−1 . The energies are plotted for different emax .
The results are independent of the underlying oscillator frequency to
∟ 1% for emax = 12. The grey dashed line is the experimental binding energy. These results underbind with respect to their doubly-magic
counterparts in Fig. 3.9. . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

Figure 3.8:

Figure 3.9:

ix

Figure 3.11: Ground-state COM energies, Eq. (3.68), for 16 O and 40 Ca at various harmonic oscillator frequencies with the NN+3N(400)-induced interaction
with λSRG = 2.0 fm−1 at emax = 12. Using the proper COM oscillator
frequencies shows the approximate factorization of Eq. (3.67). . . . . . .
Figure 4.1:

Figure 4.2:

75

The 42 lowest single-particle states (the first 5 shells) in the 2D harmonic
oscillator basis. Each box represents a single-particle state arranged by
m` , ms , and energy, and the up/down arrows indicate the spin of the
states. Within each column, the principal quantum number n increases
as one traverses upward. . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

Ground state energy (in Hartrees) of quantum dots with N particles and
an oscillator frequency of ω calculated with several different methods. .

93

Figure 4.3:

Particle-attached energy (in Hartrees) of quantum dots with N particles
and an oscillator frequency of ω calculated with several different methods. 94

Figure 4.4:

Particle-removed energy (in Hartrees) of quantum dots with N particles
and an oscillator frequency of ω calculated with several different methods. 94

Figure 4.5:

Energy difference of particle-attached and particle-removed states between the EOM-CC method and the exact FCI method for a quantum
dot with various parameters plotted against the single-particle overlap of
the FCI state, n1-particle or n1-hol , see Eqs. (4.42) and (4.43). The strong
correlation shows that the quality of EOM states can be judged by this
metric. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

Ground-state energies for the particle-attached nuclei 17 O, 17 F, 23 O,
and 23 F as a function of the harmonic oscillator energy ~ω with the
NN+3N(400) interaction, SRG softened with λSRG = 2fm−1 . The energies are plotted for emax = 8, 10, 12, showing the convergence as the
model space increases. The grey dashed line is the experimental binding
energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

Ground-state energies for the particle-removed nuclei 15 N, 15 O, 21 N,
and 21 O as a function of the harmonic oscillator energy ~ω with the
NN+3N(400)-induced interaction, SRG softened with λSRG = 2fm−1 .
The energies are plotted for emax = 8, 10, 12, showing the convergence
as the model space increases. The grey dashed line is the experimental
binding energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

98

Figure 4.6:

Figure 4.7:

x

Figure 4.8:

Ground-state COM energies, Eq. (3.68), for open-shell nuclei at varies
harmonic oscillator frequencies with the NN+3N(400)-induced interaction with λSRG = 2.0 fm−1 at emax = 12. The top row shows the results
for the particle-attached nuclei 17 O, 17 F, 23 O, and 23 F, and the bottom
row shows the results for the particle-removed nuclei 15 N, 15 O, 21 N, and
21 O. The right column shows that the COM wave function practically
vanishes according to Eq. (3.67). . . . . . . . . . . . . . . . . . . . . . . 100

Figure 4.9:

Low-lying PA-EOM states for 17 O and 17 F with and without a LawsonGloeckner term, along with the experimentally determined spectra. The
negative-partity states are COM contaminants and are removed by artificially raising the COM excitation energy with the parameter β according
to Eq. (3.72). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

Figure 4.10: Low-lying PA-EOM states for 23 O and 23 F with and without a LawsonGloeckner term. The negative-partiy states in 17 O are COM contaminants and are removed by artificially raising the COM excitation energy
with the Lawson-Gloeckner method, Eq. (3.72). The excited states of
these nuclei have not been experimentally detemined. . . . . . . . . . . 101
Figure 4.11: Low-lying PR-EOM states for 15 N and 15 O with and without a lawsongloeckner term, and the experimentally determined spectra. The 1/2+
state is a COM contaminant and is removed by artificially raising the
COM excitation energy with the parameter β according to Eq. (3.72). . 102
Figure 4.12: Low-lying PA-EOM states for 21 N and 21 O with and without a LawsonGloeckner term, Eq. (3.72), and the experimentally determined spectra.
The negative-partiy states in 17 O are COM contaminants and are removed by artificially raising the COM excitation energy. The excited
states of these nuclei have not been experimentally detemined. . . . . . 102
Figure 5.1:

Schematic representations of the three free-space weak processes in this
work: β − decay (a), β + decay (b), and electron-capture (c). The coupling
constant for the point interaction vertex is the weak-interaction coupling
constant gW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Figure 5.2:

Schematic representations of a higher-order weak interactions involving
pion-exchange that occur within a nucleus, for β − decay (a), β + decay
(b), and electron-capture (c). These processes are not included in the
impulse approximation. The coupling constant for the point interaction
vertex is the effective weak-interaction coupling constant GF . . . . . . . 105

xi

Figure 5.3:

Schematic representations of the impulse approximations to the different
weak processes within an A-body nucleus: β − decay (a), β + decay (b),
and electron-capture (c). The active nucleon doesn’t interact with the
initial and final nuclei during the weak process. The coupling constant
for each interation vertex is the weak-interaction coupling constant GF .

xii

106

KEY TO SYMBOLS AND ABBREVIATIONS

Ap -Bh – A-particle, B-hole excitation or de-excitation from the reference state
|0i – vacuum state
|Φ0 i – reference state
|Ψi – correlated ground state
a ¡¡¡a
|Φi 1···i A i – specific Ap -Bh state
1

B

{· · · } – normal-ordered with respect to the reference state
~ω – harmonic oscillator energy scale
Ĥ – Hamiltonian
H – similarity-transformed Hamiltonian
ĤN – normal-ordered Hamiltonian
HN – normal-ordered similarity-transformed Hamiltonian
T̂ – cluster operator
aA
i
a ¡¡¡a
a
i
εi 1···i A – energy denominator, fi 1 + · · · + fi B − fa11 − · · · − faA
1
1 B
B

CC – coupled cluster
CCD – coupled cluster with doubles
CCSD – coupled cluster with singles and doubles
CCSDT – coupled cluster with singles, doubles, and triples
CI – configuration interaction
FCI – full configuration interaction
COM – center of mass
EOM – equation-of-motion

xiii

PA – particle-attached
PR – particle-removed
HF – Hartree-Fock
IM-SRG – in-medium similarity renormalization group
HO – harmonic oscillator
MBPT – many-body perturbation theory
DIIS – direct-inversion of the iterative subspace
QCD – quantum chromodynamics
EFT – effective field theory
NN – nucleon-nucleon
3N – three-nucleon
NLO – next-to leading order
N2 LO – next-to-next-to leading order
N3 LO – next-to-next-to-next-to leading order

xiv

Chapter 1
Introduction
Steady progress in any modern scientific endeavor requires a strong, dynamic relationship
between experimental data to paint an accurate picture of some natural phenomena and
theoretical models to interpret those phenomena with respect to the growing network of other
scientific models. Conversely, the predictive capability of theoretical models can highlight
blurry or unfinished areas of that picture which can be clarified or completed by new or
improved experimental techniques. In the pursuit to understand and describe the atomic
nucleus and the corresponding implications from quarks to neutron stars, this push-andpull coordination between theory and experiment makes progress in modern nuclear physics
robust and persistent.
An integral component of modern nuclear physics is describing the structure and emergent properties of self-bound systems of protons and neutrons. The systems in questions can
be stable nuclei, rare isotopes far from stability, and even infinite nuclear matter which can be
used to model neutron stars. Relevant properties to nuclear structure include ground-state
energies–for determining nuclear masses, excited-state energies–for identification in gamma
or neutron spectroscopy, and transition or decay amplitudes–for calculating the respective
rates for those processes. This wide array of emergent properties inserts both nuclear structure theory and experiment into a prominent role within every other subfield of modern
nuclear physics, from lattice quantum chromodynamics (QCD) to nuclear astrophysics, and
beyond, to questions about fundamental symmetries and dark matter. However, two inextricable characteristics of a comprehensive model of nuclear structure–the increasingly large
1

size of many-body nuclear systems and the complexity and strength of the nucleon-nucleon
interactions–have been imposing hurdles for theorists to overcome.

1.1

A Brief History of Nuclear Structure Theory

A major step in the project to solve the correlation problem in many-fermion systems was
taken with the work of Brueckner, Bethe, and Goldstone [5, 6, 7] with the reformulation of the
nuclear interaction by accounting for two-body correlations from the nuclear medium. This
work continued with the work of Coester [8, 9] and Kummel [10], amongst many others, with
a further resummation of nuclear correlations in the form of an exponential ansatz into what
would become coupled-cluster (CC) theory. However, there were two major obstacles that
hindered the progress in this area for decades. First, while these methods were systematically
improvable by including progressively higher-level correlations, the highly nonperturbative
nature of the nuclear force required computationally infeasible summations. Second, there
wasn’t a reliable and consistent theory to model nucleon-nucleon interactions.
On the other hand, with the well-known and highly-perturbative Coulomb force, which
underlies the many-electron systems in atoms and molecules, the field of quantum chemistry
made consistent advances in ab-initio quantum chemistry possible since the 1950s. Along
with the quasi-exact method of configuration interaction (CI) theory [11, 12, 13, 14] which
have been utilized since the formulation of quantum mechanics, chemists successfully employed approximate methods like many-body perturbation theory (MBPT) [15, 16, 17, 18]
and coupled-cluster theory [19, 20, 21, 22, 18].
Fortunately, within the past decade, two breakthroughs have allowed ab initio nuclear
structure to resurface and thrive the way that quantum chemistry had done in the previous
decades. First was the invention of chiral effective field theory (EFT) [23, 24] which gave
2

Figure 1.1: Nuclear chart of nuclei with ground-state energies which have been calculated
with ab-initio methods and NN+3N interactions. Figure taken from [1].
theorists the ability to construct nucleon-nucleon interactions consistent with the underlying
QCD of the strong nuclear force. Second was the application of renormalization group (RG)
methods to the nuclear force [25, 26]. This procedure can “soften” the NN interaction, to
decouple the high- and low-momentum components of the nuclear force and generate lesscorrelated systems that can be calculated at a reasonable computational cost. These major
changes to nuclear structure theory made it possible to merge the field with the progress of
quantum chemistry and open a new area for additional developments in ab initio descriptions

3

of many-fermion systems, see Fig. 1.1.
Along with exponential improvements to high-performance computing, these novel techniques have allowed modern many-body methods to extend their reach and deepen their
applicability across the nuclear chart, see Fig. 1.2. The no-core shell model (NCSM), an
exact method for a given model space, has been useful in calculating the radii, transition strengths, and effective interactions of light nuclei and has been extended to nuclei
in the sd shell [27, 28, 29]. A quasi-exact technique which follows a completely different
methodology than NCSM, quantum Monte Carlo (QMC), has also progressed and is now
capable of calculating properties of light nuclei with modern chiral forces [30, 31, 32]. In
addition to these exponentially scaling techniques’ successes with lighter nuclei, polynomially scaling techniques–such as the in-medium similarity renormalization group (IMSRG)
[33, 34, 35, 36, 37, 38, 39, 40], self-consistent Green’s functions (SCGF) [41, 42, 43], and coupled cluster theory [44, 45, 46, 47, 48, 49, 50, 51]–have been able to reach open-shell nuclei
through the pf shell and even up to the chain of even tin isotopes with equations-of-motion
and multi-reference techniques [52].

1.2

Electroweak Theory and Nuclear Structure

Nuclear structure is implicated in performing and analyzing experiments to probe fundamental symmetries and physics beyond the Standard Model. One example is determining
the Vud component of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, which relates quark
eigenstates of the weak interaction to their mass eigenstates [53, 54]. This component can be
determined from by measuring the half-lives of superallowedbeta decays [55] and applying
a nucleus-dependent structure correction [56, 57, 58, 59, 60]. The value of |Vud | is used to
test the unitarity of the CKM matrix and the conserved-vector current hypothesis, which
4

140

Mass Number (A)

120
100
80
60
40
20
0

1985

1990

1995

2000

Year

2005

2010

2015

Figure 1.2: Progress of ab-initio nuclear structure from calculations of ground-state energies
with NN+3N interactions. Early progress was approximately linear as the problem size
scaled with Moore’s law while more recent progress has taken advantage of new algorithms
which have outpaced Moore’s law. Data taken from [1].
relates the f t-values of superallowed Fermi β decays of different nuclei, both predicted by
the standard model [61].
Another example of physics beyond the standard model is the neutrinoless double-beta
decay (0νββ) [62, 63]. The extremely-rare, two-neutrino double-beta decay (2νββ) has been
observed in many experiments [64, 65], which has motivated the search for its neutrinoless
counterpart, in which two Majorana neutrinos, being their own antiparticles, annihilate
one another, which is not possible in the standard electro-weak theory. The long half-lives
of these theoretical decays depend on a phase-space factor, which is highly dependent on
the decay Q-value, and a nuclear matrix element. The Q-value can be determined from
high-precision mass measurements of the relevant nuclei [66, 67, 68, 69], while the nuclear
matrix element, which contributes the largest source of uncertainty, must be calculated with

5

a reliable many-body theory.
The weak interaction and nuclear structure can also be exploited for supernova neutrino
detection and spectroscopy. While these original detectors were based on electron-neutrino
scattering [70, 71], more recent experiments utilize correlated nucleon effects of large nuclei
to enhance the scattering cross section and therefore the ability to resolve energies and
distinguish neutrino flavors [72, 73, 74, 75]. Supernova models predict distinct distributions
for different neutrino flavors based on the temperatures at which they are emitted [76, 77].
With nuclear structure calculations that include sufficient nuclear correlations, these highresolution detectors can be used to verify specific models.

1.3

Ab-Initio Descriptions of Beta Decay

Since Enrico Fermi’s originally rejected paper describing beta decay in 1934 [78, 79], theorists have worked to refine this description within the ever-growing library of knowledge
concerning the nature of the weak force, the characteristics of the neutrino, and the structure
of nuclei. With the success of ab initio calculations for nuclear properties such as masses,
radii, and electromagnetic phenomena, these techniques also seem promising ways to calculate relevant quantities involved in nuclear beta decay. Because the kinematics of the decay
and the underlying weak process are well understood, the remaining task for nuclear theory
to tackle is calculating the transition amplitudes between the initial and final nuclei.
Modern calculations of these beta-decay matrix elements were originally performed using
phenomenological interactions in the shell model framework [80, 81, 82, 83]. Also, predecessors to current ab-initio techniques like the random-phase approximation (RPA) [84] included
core-correlation effects in these early descriptions. These methods were able to successfully
reproduce experimental lifetime data and address technical issues such as the quenching of
6

the axial-vector coupling constant. More recently, the success of the shell model has inspired
an extension to the new method, known as the ab-initio shell model [29], where an effective
interaction is constructed within a certain valence space using a many-body method such as
CC [85] or IMSRG [86]. However, these techniques are computationally expensive and cannot currently reach heavy nuclei of interest. The most common method used in their place
is known as the quasiparticle random-phase approximation (QRPA) [87, 88]. While these
calculations can be performed for heavy nuclei in large spaces, they also rely on phenomenological effective interactions. Therefore, there is a demand for computationally-economical,
ab initio techniques that can capture the relevant many-body correlations needed to accurately describe the nuclear structure aspects of electro-weak processes.

1.4

Thesis Structure

The main goal of this work is to explore the ab initio description of nuclear beta decay within
the coupled-cluster theory framework of the equation-of-motion coupled cluster with singles
and doulbes (EOM-CCSD) method using renormalized chiral NN and 3N interactions. The
organization of the thesis builds from a general description of the many-body problem of
quantum mechanics in chapter 2. Then, in chapter 3, this many-body framework is applied
within the coupled-cluster theory and applied to various systems including, atomic nuclei. In
chapter 4, coupled-cluster theory is extended to the equation-of-motion method to describe
open-shell systems. Then, in chapter 5, different types of beta decay are described in detail
then outlines the procedure to express beta-decay observables as effective coupled-cluster operators and how to calculate those observables in the equation-of-motion framework. Lastly,
conclusions and future perspectives are given in chapter 6 while technical details concerning
the formalism and implementation are given in the appendices A – D.
7

Chapter 2
Many-Body Quantum Mechanics
Ab initio structure calculations of many-fermion systems such as those in nuclear and electronic structure aim to describe emergent phenomena from the constituent particles subject
to the underlying microscopic Hamiltonian. This amounts to finding the solution to the
many-body Schrödinger equation. However, a calculation of the exact solution needs to
account for all possible correlations among the particles and thus scales factorially. This motivates the need for approximations to the exact solution that account for the most important
correlations. This chapter first establishes the formalism necessary to define the many-body
problem then illustrates several successive approximations to its solution. Because the type
of fermions and the underlying Hamiltonian can be kept generic until specific systems are
considered, the formalism and many-body methods can be kept generic as well.

2.1

Independent-Particle Model

The nonrelativistic A-body quantum problem begins with the Schrödinger equation,

ĤΨν (r1 , · · · , rA ) = Eν Ψν (r1 , · · · , rA ) ,

(2.1)

for the correlated wave function Ψν (r1 , ¡ ¡ ¡ , rA ) and the corresponding energy Eν . The
Hamiltonian can be written generically as a sum of k-body pieces which, in principle, can

8

contain up to A-body interactions,

Ĥ = (1)Ĥ + (2)Ĥ + (3)Ĥ + · · ·
=

A
X
i

(1)Ĥ(r ) +
i

A
X

(2)Ĥ




ri , rj +

i<j

A
X

(3)Ĥ

i<j<k



ri , rj , rk + ¡ ¡ ¡ .

(2.2)

2

2
The one-body term can contain the kinetic energy operator, −~
2m ∇i , as well as any external

potential while the higher-order terms contain inter-particle interactions.
An intuitive way to formulate the solution to the many-body Schrödinger equation is to
express the collective wave function in terms of independent single-particle wave functions, or
orbitals φ (r). In this independent-particle model, a selection of single-particle wave functions,
known as the single-particle basis, are constructed by solving the Schrödinger equation for a
single particle in either a mean-field potential for bound systems or in free space for infinite
systems. Then a many-body wave function is constructed as a product of these singleparticle orbits. This simple model is justified because it becomes exact when inter-particle
interactions are completely suppressed and is useful because it provides an intuitive way to
interpret complicated many-body dynamics as processes involving few single-particle wave
functions.
A many-body wave function of fermions must be anti-symmetric with respect to particle
exchange so that the Pauli exclusion principle is followed, such that no single-particle wave
function is occupied by more than one fermion. This condition is satisfied by a wave function
in the form of a Slater determinant [89],

 φ1 (r1 ) φ1 (r2 )

1  φ2 (r1 ) φ2 (r2 )
Φ (r1 , · · · , rA ) = √  .
..
A!  ..
.
φ (r ) φ (r )
A 1
A 2
9

¡¡¡
¡¡¡
...
¡¡¡


φ1 (rA ) 
φ2 (rA ) 
,
..

.

φ (r )
A

A

(2.3)



where A is the number of particles in the system and φp rµ is the p-th orbital filled with
the Âľ-th particle.
If the orbitals are constructed from an appropriate phenomenological potential, a Slater
determinant composed of the A lowest orbitals can represent a fairly good approximation
to the ground state for a closed-shell system, where the lowest-energy Slater determinant
can be uniquely determined. The set of all Slater determinants in a certain model space of
single-particle wave functions defines a complete A-body Hilbert space such that a generic
wave function can be written as a linear combination of Slater determinants,
N
X
Âľ
Ψν (r1 , ¡ ¡ ¡ , rA ) =
Cν Ό¾ (r1 , ¡ ¡ ¡ , rA ) ,

(2.4)

Âľ=1

Âľ

Âľ

where Cν = hΨ (r1 , ¡ ¡ ¡ , rA ) |Όν (r1 , ¡ ¡ ¡ , rA )i. The number of Slater determinants N in an
A-body Hilbert space with N orbits is given by,

N =

N
A

!
=

N!
,
A!(N − A)!

(2.5)

which shows the factorial scaling of the exact problem. However, to reduce the size of the
problem, progressively more significant Slater determinants can be chosen to systematically
refine approximations to the full solution.

2.2

Second Quantization

Even with the simplification of the independent-particle model, the many-body Schrödinger
equation is an unwieldy and complex system of coupled differential equations. A useful
reformulation of this equation is to promote the single-particle orbits to operators in a step

10

known as second quantization (see e.g., [18, 90]). In this framework, a Slater determinant is
represented by a string of occupied orbitals,



Φ (r1 , · · · , rA ) ≡ A φp1 φp2 φp3 · · · φpN



≡ |p1 p2 p3 · · · pN i,

(2.6)

where A represents a permutation and normalization operator to correspond with Eq. (2.3).
These second-quantized Slater determinants can be constructed with the use of operators
†

that correspond to specific orbitals. A creation operator, âp , places a particle in the p orbital,
and an annihilation operator, âp , removes a particle from the p orbital,

†

âp |0i = |pi

âp |pi = |0i,

(2.7)

where |0i represents the true vacuum, a state void of any particles. Because there must
be a correspondence between the original first quantization and second quantization, these
creation an annihilation operators obey the following anticommutation relations ([Â, B̂]+ =
ÂB̂ + B̂ Â),
†

[âp , âq ]+ = δpq

†

†

[âp , âq ]+ = [âp , âq ]+ = 0,

(2.8)

which guarantee that wave functions comprised of these operators obey antisymmetry and
the Pauli exclusion principle required of fermionic systems.
The Hamiltonian in the form of Eq. (2.2) can be written with second-quantized operators
as,

Ĥ =

X
pq

(1)H p
q

†

âp âq +

1 X (3) pqr † † †
1 X (2) pq † †
Hrs âp âq âs âr +
Hstu âp âq âr âu ât âs + · · · , (2.9)
4 pqrs
36
pqrstu

where the prefactors account for the double counting of particle-particle interactions, and

11

the matrix elements represent integrals over the relevant single-particle wave functions,

(1)H p
q
(2)H pq
rs

≡

Z

≡

Z

dr1 φ∗p (r1 ) (1)Ĥ(r1 ) φq (r1 )

dr1 dr2 φ∗p (r1 ) φ∗q (r2 ) (2)Ĥ(r1 , r2 ) [φr (r1 ) φs (r2 ) − φs (r1 ) φr (r2 )]
..
.

(2.10)

Matrix elements involving two or more particles include exchange terms which guarantee
that they are also antisymmetric,

(2)H pq
rs
(3)H pqr
stu

qp

pq

qp

= −(2)Hrs = −(2)Hsr = (2)Hsr
pqr

qpr

qpr

= −(3)Hstu = −(3)Htsu = (3)Htsu = · · ·

(2.11)

These definitions apply regardless of the form of the Hamiltonian, and thus this formalism
remains generic to the particular system. Second quantization is a crucial step in simplifying the many-body Schrödinger equation because it reduces the complexity of the spatial
and spin degrees of freedom within the single-particle wave functions and interactions into
precomputed matrix elements. The remaining effort is reduced to algebraic expressions involving creation and annihilation operators.

2.3

Normal Ordering

It’s convenient to define an A-particle reference state, where states are filled from the true
vacuum up to a closed shell, known as the Fermi level. This reference state must be uniquely
determined from the number of particles in the system and therefore nondegenerate with

12

|ÎŚ0 i =

EF

Figure 2.1: A depiction of the closed-shell reference state in the independent particle model.
Each horizontal line represents a shell of single-particle orbits, represented by circles, and the
dotted line represents the Fermi level which separates the unoccupied particle states from
the occupied hole states.
other Slater determinants,
|ÎŚ0 i =


A
Y


i

†

âi





|0i.

(2.12)

This reference determinant defines a new Fermi vacuum. States above the Fermi vacuum are
called particle states and will be denoted with the indices a, b, c, d... while states below the
Fermi vacuum are called hole states and will be denoted with the indices i, j, k, l.... Generic
states above or below the Fermi vacuum will be denoted with the indices p, q, r, s....
Any other Slater determinant can be constructed relative to this reference state by adding
particles and/or removing holes. A Slater determinant with A particles added and B holes
removed from reference state is known as a Ap -Bh excitation. A 1p -1h state is constructed
by removing a particle in the occupied state i and adding a particle in the unoccupied state
a,
†

|Φai i ≡ âa âi |Φi.

(2.13)

Equivalently, a 2p -2h state is constructed by removing particles in states i and j then adding
them to states b and a,
† †

|ÎŚab
ij i ≡ âa âb âj âi |Φi.

13

(2.14)

|ÎŚai i =

|ÎŚab
ij i =

EF

|ÎŚa i =

|ÎŚi i =

EF

EF

EF

Figure 2.2: A depiction of 1p -1h, 2p -2h, 1p -0h, and 0p -1h Slater determinants defined
relative to the reference state in the independent particle model.
The number of creation and annihilation operators doesn’t neccessarily have to be equal. For
instance, a single particle can be added on top of the reference state with a single creation
operator,
†

|Φa i ≡ âa |Φi,

(2.15)

and a particle can be removed with a single annihilation operator,

|Φi i ≡ âi |Φi.

(2.16)

Using these definitions, hole-creation and particle-annihilation operators vanish when act†

ing on the Fermi vacuum from the left, âi |Φ0 i = âa |Φ0 i = 0. Conversely, hole-annihilation
and particle-creation operators vanish when acting on the Fermi vacuum from the right,
†

hΦ0 | âi = hΦ0 | âa = 0.
These results can be exploited to simplify expressions involving strings of creation and
annihilation operators by a procedure called normal ordering with respect to the Fermi
vacuum. Denoted by {¡ ¡ ¡ }, normal ordering permutes a string of creation and annihilation
14

operators so that hole-annihilation and particle-creation operators are to the left of holecreation and particle-annihilation operators, which guarantees that normal ordered operators
vanish on the Fermi vacuum, hΌ0 | {¡ ¡ ¡ } = 0 and {¡ ¡ ¡ } |Ό0 i = 0.
n
o
†
†
†
†
âj · · · âi · · · âb · · · âa = (−1)σ âi · · · âa · · · âj · · · âb ,

(2.17)

where σ is the number of two-state permutations required to do the normal ordering.

2.4

Wick’s Theorem

At this point, the many-body problem has been reduced to computing long strings of creation
and annihilation operators between the normal-ordered Hamiltonian and the correlated wave
function using Eq. (2.8). Instead of using a brute-force approach by permuting over and over,
a further simplification known as Wick’s theorem [91] can be introduced. A Wick contraction
of two operators with respect to the reference state is defined as

n
o
ÂB̂ = ÂB̂ − ÂB̂ .

(2.18)

Which, given the definition in Eq. (2.17) and the anticommutation relations int Eq. (2.8),
means that the only nonzero contractions are of the form,

†

âi âj = δij

and

†

âa âb = δab .

(2.19)

Because contracted operators simply represent a Kronecker delta, they can be removed from
a normal ordered product by permuting the product σ times so that the contracted operators

15

are next to each other,

{Â · · · B̂ · · · Ĉ · · · D̂} = (−1)σ {Â · · · B̂ Ĉ · · · D̂} = (−1)σ B̂ Ĉ{Â · · · D̂}.

(2.20)

These different definitions for operator manipulation come together to define the timeindependent Wick’s theorem, which reformulates a product of operators as the sum of its
normal-ordered form and all possible contractions of its normal-ordered form,

ÂB̂ Ĉ · · · =

n
o
X
X
X
{ÂB̂ Ĉ · · · } +
{ÂB̂ Ĉ · · · } + · · · +
{ÂB̂ Ĉ · · · }.
ÂB̂ Ĉ · · · +
onecontractions

twocontractions

allcontractions

(2.21)
Wick’s theorem is incredibly useful in many-body techniques because complicated expressions of operators can be expressed as diagrams that are easy to compute with simple
diagrammatic rules which correspond to Eqs. (2.8),(2.19), and (2.20). These diagrammatic
techniques are an integral component to deriving expressions used in this work, and their
underlying rules are extensively discussed in [18].
A powerful application of Wick’s theorem is to rewrite the Hamiltonian in Eq. (2.2) as a
sum of normal-ordered operators in the form of Eq. (2.21),

Ĥ = E0 +

o
o
X p n † o 1 X pq n † †
1 X pqr n † † †
fq âp âq +
Vrs âp âq âs âr +
Wstu âp âq âr âu ât âs + · · · .
4 pqrs
36
pq
pqrstu

(2.22)
This form of the Hamiltonian can be split into a zero-body component E0 , known as the
reference energy, and the remaining normal-ordered Hamiltonian, ĤN ,

ĤN =

o
o
X p n † o 1 X pq n † †
1 X pqr n † † †
Vrs âp âq âs âr +
Wstu âp âq âr âu ât âs + · · · .
fq âp âq +
4 pqrs
36
pq
pqrstu

(2.23)

16

The reference energy, E0 , contains fully-contracted terms, and because the Hamiltonian operators are ordered so that the creation operators appear before the annihilations
†

operators, only terms that contract hole states in the form {· · · âi · · · âj · · · } are nonzero.
Therefore, the zero-body component of the normal ordered Hamiltonian can be written as a
sums over all hole states for each component of the original Hamiltonian,

E0 =

X

(1)H i
i

+

1 X (2) ij 1 X (3) ijk
Hij +
Hijk ¡ ¡ ¡ .
2
6
ij

i

(2.24)

ijk

In a compact diagrammatic form, this sum can be drawn as the sum of connected vertices
(

,

, etc...) corresponding to the components of the original Hamiltonian in

Eq. (2.2),
i

E0 =

i

j

+

i

j

k

+ ¡¡¡ .

+

(2.25)

The number of lines connected to each vertex defines its type such that the single vertex
corresponds to the one-body Hamiltonian, the double vertex corresponds to the two-body
Hamiltonian, etc... The looped lines represent a sum over all hole states. The reference
energy is also equivalent to the Hamiltonian expectation value for the reference state,

E0 = hΦ0 | Ĥ|Φ0 i.

(2.26)

p

The second component of the normal-ordered Hamiltonian, fq

n
o
†
âp âq , contains terms

that have all operators contracted except for one pair. Like the fully-contracted term, these

17

contractions must involve only hole states,

p

p

fq = (1)Hq +

X

(2)H pi
qi

1 X (3) pij
Hqij + ¡ ¡ ¡ .
2

+

(2.27)

ij

i

In diagrammatic notation, the uncontracted pair of operators is represented as two external
lines connected to each vertex which, because they are generic states, are drawn as unoriented
lines,

p

p

p

=
q

p

i

+

j

i

+ ¡¡¡ .

+

q

q

q

(2.28)

When a line represents a particle state, it will contain an arrow directed upwards while a
line representing a hole state will contain an arrow directed downwards.
n
o
pq
† †
The two-body term, Vrs âp âq âs âr , follows in the same manner such that it represents
all terms of the original Hamiltonian with all but two pairs of operators contracted. Like the
zero-body and one-body terms, the two-body term contains the original two-body term (2)Ĥ
as well as density-dependent terms that sum over hole states in higher-body Hamiltonian
terms, leaving four external lines for each diagram vertex,

pq

pq

Vrs = (2)Hrs +

X

(3)H pqi
rsi

i
p

q

p

q

=
r

s

+ ¡¡¡
p

q

+
r

s

r

s

i

+ ¡¡¡ .

(2.29)

Three- and four-body normal-ordered terms can also be calculated by following the same

18

procedure, but in practice are truncated at the two- or three-body level. Because normal
ordering the Hamiltonian has the effect of shuffling higher-order interactions into lower-order
terms, it becomes feasible to include computationally expensive many-body interactions as
normal-ordered few-body interactions. Also, it reorganizes many-body correlations into the
reference state so that additional correlations around the Fermi surface can be treated as a
perturbation. Therefore, from this point forward, the many-body problem will be formulated
in terms of the normal-ordered terms in Eq. (2.22), and the bare interactions will be truncated
beyond the three-body level for computational feasibility. Electronic systems are naturally
truncated at the level of the two-body Coulomb interaction, while nuclear systems can be
successfully described with the two-body normal-ordered piece of the three-body force, in
the form of Eq. (2.29).
With this new partition, the many-body Schrödinger equation for the ground state |Ψi
can be written in terms of the normal-ordered Hamiltonian as,

Ĥ|Ψi = (E0 + ĤN )|Ψi = E |Ψi
−→ ĤN |Ψi = (E − E0 )|Ψi = ∆E |Ψi,

(2.30)

where ∆E is known as the correlation energy.
Now that the many-body quantum problem has been formulated, different approaches to
solving that problem can be proposed and analyzed. Because taking account of correlations
from all particles simultaneously is a demanding–and for some systems, computationally
impossible–endeavor, methods for solving the many-body Schrödinger equation should be
systematically improvable. Successful methods with this quality incorporate the most dominant correlations in lower-order solutions and approach the exact solution when more and

19

more orders are included.

2.5

Hartree–Fock Method

A successful, first-order approximation to any many-body method comes from noticing that
each individual particle feels a mean-field potential from the cumulative interactions with all
the other particles. The Hartree-Fock (HF) method [92, 93] aims to transform the original
single-particle basis to a Hartree-Fock basis where each orbital is the eigenfunction of its
corresponding mean-field. Because the transformation of a single orbital changes its effect on
every other particle, this process must be performed iteratively until self-consistency between
all the orbitals is reached, which is why this method is also known as the Self-Consistent
Field (SCF) method.
This mean-field picture results from the following procedure. It begins by minimizing the
reference energy with respect to the reference state. This functional is just the zero-body
piece of the normal-ordered Hamiltonian,

EHF [Φ0 ] = hΦ0 | Ĥ|Φ0 i =

X

(1)H i
i

i

+

1 X (2) ij 1 X (3) ijk
Hij +
Hijk .
2
6
ij

(2.31)

ijk

Transforming the reference determinant can be accomplished by rotating the state within the
single-particle basis by use of the Thouless theorem [94], which states that any Slater determinant can be written as the product of any other Slater determinant and an exponentiated
single-excitation operator,

|Φ0 i = eĈ1 |Φ0 i,

where Ĉ1 =

X
ai

20

n
o
†
Cia âa âi .

(2.32)

If the difference between the two Slater determinants is dominated by single excitations, this
transformation can be approximated by expanding the exponential and ignoring higher-order
terms,
|ÎŚ0 i

'

1+

X

Cia

ai

!
n
o
†
âa âi
|ÎŚ0 i.

(2.33)

The reference energy functional can now be written as a sum of the original reference
state and new terms that incorporate the single-excitation variation,

X
X
 
a
Cia hΦ0 | Ĥ|Φai i +
C a∗
EHF Φ0 = hΦ0 | Ĥ|Φ0 i ' EHF [Φ0 ] +
i hΦi | Ĥ|Φ0 i.
ai

(2.34)

ai

The minimum of this functional is found by differentiating the expression with respect to
the coefficients Cia and setting the result to zero,

X
  X a
a
δC a∗
δEHF Ό0 '
δCi hΦ0 | Ĥ|Φai i +
i hΦi | Ĥ|Φ0 i = 0.

(2.35)

ai

ai

Because this expression is Hermitian, both terms must vanish independently so that,

hΦ0 | Ĥ|Φai i = hΦai | Ĥ|Φ0 i = 0.

(2.36)

This condition is the result of the Brillouin theorem [95], which states that the Hamiltonian
matrix element must vanish between an optimized Hartree-Fock ground state and any single
excitation from it. The Brillouin condition is satisfied by diagonalizing the one-body piece
p

of the normal-ordered Hamiltonian fq (Eq. (2.27)), known as the Fock operator, such that
off-diagonal pieces like hΦ0 | Ĥ|Φai i = fai and hΦai | Ĥ|Φ0 i = fia vanish. Diagonalizing the Fock

21

operator can be written schematically as,

p

ξq δpq

p

q

p
←− (1)Hq +

p

←−

X

(2)H pi
qi

1 X (3) pij
Hqij
2
ij

i
p

p

i

+
q

+

i

j

+
q

,

(2.37)

q

p

where Îľq is the eigenvalue of the Fock operator, and its diagrammatic form vanishes when
the external indices differ.
A practical way of solving this system of equations is to express each new orbital in
the unknown Hatree-Fock basis, |p0 i ≡ φp0 (r), denoted with primed labels, as a linear
combination of the known single-particle basis states, |pi ≡ φp (r), denoted without primed
p

labels. These two bases are related by a unitary transformation C 0 ≡ hp|p0 i,
p
|p0 i =

X
p

hp|p0 i|pi =

X p
C 0 |pi.
p

p

(2.38)

Then the Fock matrix can be written in terms of the Hartree-Fock basis,
p0

p0

f 0 = (1)H 0 +
q
q

X
i0

(2)H p0 i0
q 0 i0

+

1 X (3) p0 i0 j 0
H 00 0
qij
2 0 0
ij

X p0 ∗
X p0 ∗ 0
1 X p0 ∗ i0 ∗ j 0 ∗ (3) prs q t u
pr q
p q
=
Cp (1)Hq C 0 +
Cp Cri ∗ (2)Hqs C 0 Cis0 +
Cp Cr Cs
Hqtu C 0 Ci0 Cj 0 .
q
q
q
2
pq
0
0 0
i
prqs

ij
prsqtu

(2.39)
p

Defining the first-order density matrix Îłq as the product of expansion coefficients, summed

22

over all shared hole states,
p

Îłq =

X p 0
C 0 Cqi ∗ ,

(2.40)

i

i0

Eq. (2.38) is simplified to,
p0
f 0
q

#
"
X
X p0 ∗
X
1
prs
q
pr
(1)H p +
Îłtr Îłus (3)Hqtu C 0
=
Cp
Îłsr (2)Hqs +
q
q
2
pq
rs
rstu

p0

−→ ε 0 δp0 q0 .
q

(2.41)

Therefore, the Hartree-Fock equations are ultimately expressed as an eigenvalue problem
according to Eq. (2.5) where the matrix to diagonalize is the Fock matrix in the form,

 
X
1 X t u (3) prs
p
p
pr
Îłr Îłs Hqtu ,
F̂q Ĉ = (1)Hq +
Îłrs (2)Hqs +
2
rs

(2.42)

rstu

p

and the matrix of coefficients, Ĉ = C 0 , is the unitary operator that transforms the matrix
p
to a diagonal form,
X p0 ∗ p   q
p0
Cp Fq Ĉ C 0 = ε 0 δp0 q0
q

pq

q

(2.43)

The iterative nature of the solution comes from the dependence of the Fock matrix on the
transformation coefficients. These Hartree-Fock equations are solved numerically by using
an iterative algorithm where the Fock matrix is built using a known set of coefficients and
diagonalized to obtain an updated set of coefficients. This process is repeated until the
unitary set of coefficients is unchanged within a certain tolerance. For most calculations,
using the identity matrix as an initial guess for the coefficients is sufficient. To improve
the rate of convergence, techniques such as the direct inversion of the iterative subspace
(DIIS) [96, 97] or Broyden’s method [98] can be implemented. And, to avoid any oscillatory
behavior around the solution, techniques such as the level-shifting method or ad hoc linear
mixing can be implemented to dampen the large changes between iterations.

23

To make use of the HF solution as the reference state for post-HF calculations, the
Hamiltonian matrix elements must be transformed to the new basis and the normal-ordered
pieces redefined to account for the additional reordering of one-particle correlations into the
HF energy. The one-body piece from Eq. (2.27) is simply the resulting eigenvalues of the
diagonalized Fock matrix,
p0

p0

f 0 = ξ 0 δp0 q0 .
q
q

(2.44)

For the two-body term, Eq. (2.29), first the density-dependent component of the three-body
interaction is written with the first-order density matrix. Then the remaining states are
transformed according to Eq. (2.38),
p0 q 0
V0 0
rs


X p0 ∗ q0 ∗ 
pq
pqt u
(2)
(3)
Hrs + Hrsu Îłt Crr0 Css0 .
=
Cp Cq

(2.45)

pqrs

The reference energy from Eqns. (2.24) and (2.26) can be written in terms of the transformed
one- and two-body Hamiltonian, Eqns. (2.27) and (2.29), and the original three-body term
using first-order density matrices,

E0 =

X 0 1 X i0 j 0 1 X
(3)H pqr Îł s Îł t Îł r .
εii0 −
V0 0+
stu p q u
i
j
2 0 0
6
0
i

(2.46)

pqrstu

ij

Additionally, any operators that are constructed in the original basis must be transformed
in a similar manner. For example, a one-body operator Ô in the Hartree-Fock basis is,

Ô =

X
p0 q 0

p0

O 0
q

n
o X 0
n
o
†
p∗ p q
†
â 0 âq0 =
Cp Oq C 0 â 0 âq0

−→

p

q

p0 q 0 pq

p0

O 0 =
q

X p0 ∗ p q
Cp Oq C 0 .
pq

24

q

p

(2.47)

Because the Hartree-Fock basis is diagonal in the one-body piece of the Hamiltonian,
any terms that include off-diagonal elements automatically vanish, greatly simplifying any
post-Hartree-Fock methods. From this point on, all calculations will use the Hartree-Fock
basis unless stated otherwise, and prime symbols will be omitted.

2.6

Configuration-Interaction

The most generic way to write a correlated wave function in a given basis is as a linear
combination of all possible Slater determinants. This expansion can, in principle, consist of
the 0p -0h reference state and all possible N p -N h excitations up to Ap -Ah excitations,
N
X


A 
X
1 2 X
a ...a
a ...a
Ci 1...i N |ÎŚi 1...i N i.
|Ψν i =
Cνi |Όνi i = C0 |Ό0 i +
1 N
1 N
N ! a ...a
νi
N =1
1 N

(2.48)

i1 ...iN

Using this form of the wave function in Eq. (2.48), the normal-ordered Schrödinger equation
can be reformulated as a standard matrix eigenvalue problem,

ĤN |Ψν i = ∆Eν |Ψν i −→ hΨµ | ĤN |Ψν i = ∆Eν hΨµ |Ψν i
=

X
¾i νi

Cµ∗ hΦµi | ĤN |Φνi iCνi = ∆Eν
i

X
¾i νi

Cµ∗ Cνi δµi νi
i



−→ CT
hÎŚ
|
Ĥ
|ÎŚ
i
−
∆E
I
Cν = 0.
Âľ
ν
ν
N
Âľ
i
i

(2.49)

In this case, the matrix elements are Hamiltonian terms that connect two Slater determinants, and the eigenvectors are the ground and excited states in the form of Eq. (2.48).
The matrix elements can be found with the help of the Slater-Condon rules [89, 99] which,
because the Hamiltonian is restricted to one- and two-body terms, require that any terms

25

connecting Slater determinants which differ by more than two single-particle states vanish.
Also, because the one-body Hamiltonian is diagonal in the Hartree-Fock basis, it only contributes to diagonal elements of the CI matrix. Some examples of these matrix elements
are,

hΦai | Ĥ|Φai i = εa − εi − Via
ia ,
ab
cd
hÎŚab
ij | Ĥ|Φij i = Vcd ,
abd
lc
hÎŚabc
ijk | Ĥ|Φijl i = −Vkd .

(2.50)

Because the configuration-interaction method exhaustively captures all the correlations
of a many-body system, it is considered an “exact” method within a certain model space and
becomes truly exact as the number of single-particle states is increased to infinity. However,
there is a price to pay for this exactness. The number of Slater determinants in a certain
model space, N , scales factorially according to Eq. (2.5) and the configuration-interaction
matrix scales as N 2 . For sufficiently-sized model spaces, the memory required for this matrix
quickly becomes unmanagable even for the largest supercomputers, see Fig. 2.3.
However, for a reference state that is a good approximation to the true ground state,
few-body excitations generally dominate the wave functions for low-lying states [100]. This
can be exploited by truncating the expansion in Eq. (2.48). Owing to the two-body nature
of the interaction, the lowest appropriate truncation is also at the two-body level, known as
configuration interaction with singles and doubles (CISD),

|Ψν i = C0 |Ό0 i +

X
ai

Cia |ÎŚai i +

1 X ab ab
Cij |ÎŚij i.
4

(2.51)

abij

This is a very straightforward and tractable way to approximate the many-body Schrödinger
26

Figure 2.3: Scaling of the matrix size and number of non-zero matrix elements for nuclear CI
calculations of light nuclei. Even for modestly-sized model spaces, the memory requirements
approach the limit of petascale supercomputers (∟ 1010 ). Figure taken from [2].
equation, and it can be systematically improved by adding more excitations such as triples
(CISDT) or triples and quadruples (CISDTQ). But the drawback to this simplicity is that
any truncated CI method is not size-extensive such that any extensive property of a system,
like the energy, would scale with the size of the system. A desirable many-body method
will be both systematically improvable and size-extensive while maintaining computational
feasibility.

2.7

Many-Body Perturbation Theory

One many-body method that is both size-extensive and systematically improvable treats
particle-particle interactions as a perturbation to the mean-field potential and is known as
many-body perturbation theory (MBPT) [101, 15, 16, 18]. The Hamiltonian is partitioned

27

into a diagonal piece and the interaction piece,

Ĥ = Ĥ0 + V̂ , with
X pn † o
Ĥ0 = E0 +
fp âp âp
and
p

o
1 X pq n † †
Vrs âp âq âs âr .
V̂ =
4 pqrs

(2.52)

When not in the Hartree-Fock basis, the interaction piece has the additional off-diagonal
n
o
P
p
†
Fock term, p6=q fq âp âq . This means that the reference state is an eigenstate of the
zero-order piece of the Hamiltonian,

Ĥ0 |Φ0 i =

E0 +

X
i

!
n
o
†
fii âi âi
|ÎŚ0 i =

!
E0 +

X

Îľi

i

(0)

|ÎŚ0 i = E0 |ÎŚ0 i.

(2.53)

Using intermediate normalization, which sets hΌ0 |Ψi = 1, the Schrödinger equation, Eq.
(2.1) for the ground state becomes,

hΦ0 | (Ĥ0 + V̂ )|Ψi = hΦ0 | Ĥ0 |Ψi + hΦ0 | V̂ |Ψi = EhΦ|Ψi
= E (0) hΦ|Ψi + hΦ0 | V̂ |Ψi = E (0) + ∆E0 = E,

(2.54)

where the energy difference is ∆E0 ≡ hΦ0 | V̂ |Ψi.
Next, the projection operators P̂ and Q̂ can be introduced,

P̂ = |Φ0 ihΦ0 |,
Q̂ =

X
n6=0

|Φn ihΦn | = 1 − |Φ0 ihΦ0 |.

(2.55)
(2.56)

The P̂ operator isolates the reference-state component of any Slater determinant while the
28

Q̂ operator isolates all components except the reference-state component out of any Slater
determinant. Both these operators are idempotent, which means that P̂ 2 = P̂ and Q̂2 = Q̂,
and because of intermediate normalization, the correlated wave function can be written as
|Ψi = (P̂ + Q̂)|Ψi = |Φi + Q̂|Ψi. Also, both operators commute with the unperturbed
part of the Hamiltonian, Ĥ0 P̂ = P̂ Ĥ0 and Ĥ0 Q̂ = Q̂Ĥ0 . These identities can be applied
to an alternate version of the Schrödinger equation which defines a particular version of
perturbation theory known as Raleigh-Schrödinger perturbation theory (RSPT) [102, 103].
In this version, the zeroth-order energy E (0) is added to both sides of the Schrödinger
equation. Acting with Q̂ and rearranging terms gives,

Q̂ (E (0) − Ĥ0 )|Ψi = Q̂ (E (0) + V̂ − E)|Ψi
Q̂ (E (0) − Ĥ0 ) Q̂|Ψi = Q̂ (V̂ − ∆E0 )|Ψi,

(2.57)

where ∆E0 ≡ E − E (0) = hΦ0 | V̂ |Ψi. The operator Q̂ (E (0) − Ĥ0 ) Q̂ is invertible because
(E (0) − Ĥ0 )−1 is never singular in Q-space. Therefore, the operator R̂0 = Q̂ (E (0) − Ĥ0 )−1 Q̂,
known as the resolvent, can be applied to both sides which gives,

Q̂|Ψi = R̂0 (V̂ − ∆E0 )|Ψi.

(2.58)

The left-hand side of this equation can be rewritten as Q̂|Ψi = |Ψi − |Φ0 i to result in the
generating equation for RSPT,

|Ψi = |Φi + R̂0 (V̂ − ∆E0 )|Ψi.

(2.59)

Because the single-particle states are eigenfunctions of the zeroth-order Hamiltonian, they

29

are also eigenfunctions of the resolvent, with the resulting eigenvalues, Îľ, are known as energy
denominators. Applying the resolvent operator to any state orthogonal to the reference state,
see Eqs. (2.13) - (2.16), gives the following relation,

1
a ¡¡¡a
a ¡¡¡a
R̂0 |Φi 1···i N i = a ···a |Φi 1···i N i, where
1 N
1 N
ξi 1¡¡¡i N
1

N

a ¡¡¡a
εi 1···i N = εi1 + · · · + εiN − εa1 − · · · − εaN .
1

N

(2.60)

Equation (2.59) can be iterated infinitely to give the solution for the fully correlated wave
function,
|Ψi =

∞ h
X
n=0

R̂0 (V̂ − ∆E0 )

in

|ÎŚ0 i.

(2.61)

Applying this form of the correlated wave function into Eq. (2.54) results in the energy
difference,
∆E0 = hΦ0 | V̂ |Ψi =

∞
X

in
hΦ0 | V̂ R̂0 (V̂ − ∆E0 ) |Φ0 i
h

(2.62)

n=0

The immediate problem with these equations is that the right-hand sides contain the
target energy difference ∆E0 for which these equations are meant to solve. This can be
remedied by expanding the right-hand sides and rearranging terms. Using the fact that
R̂0 ∆E0 |Φ0 i = ∆E0 R̂0 |Φ0 i = 0, the first-order energy E (1) = hΦ0 | V̂ |Φ0 i, and the shifted
term Ve ≡ V̂ − E (1) , these simplify to,
|Ψi − |Φ0 i = R̂0 V̂ |Φ0 i + R̂0 Ve R̂0 V̂ |Φ0 i
+ R̂0 Ve R̂0 Ve R̂0 V̂ |Φ0 i − hΦ0 | V̂ R̂0 V̂ |Φ0 iR̂02 V̂ |Φ0 i + · · ·

(2.63)

∆E0 = hΦ0 | V̂ |Φ0 i + hΦ0 | V̂ R̂0 V̂ |Φ0 i + hΦ0 | V̂ R̂0 Ve R̂0 V̂ |Φ0 i
+ hΦ0 | V̂ R̂0 Ve R̂0 Ve R̂0 V̂ |Φ0 i − hΦ0 | V̂ R̂0 V̂ |Φ0 ihΦ0 | V̂ R̂02 V̂ |Φ0 i + · · ·
30

(2.64)

The order of each term can be easily identified by counting the numbers of times that V̂
or Ve appears. At the third order in the wave function and the fourth order in the energy,
renormalization terms make their first appearance. These terms contain separated and closed
factors in the form of lower-order energy terms, such as hΦ0 | V̂ R̂0 V̂ |Φ0 i ≡ E (2) . Terms that
do not contain normalization factors are known as principal terms.

2.7.1

Factorization Theorem

A powerful application of diagrammatic techniques known as the factorization theorem [16,
104, 105] can immediately be used to simplify these expansions. By factoring sums of unlinked
diagrams, where two or more parts of a diagram are closed and separated, from the principal
terms, it can be shown that they exactly cancel with the renormalization terms at each
order. In the following factorization, two fourth-order energy diagrams which differ by only
the time-ordering of the interaction vertices are added together. By multiplying each term
by an appropriate factor so that they share a common denominator, the additive property
cd
abcd
of the energy denominators (Îľab
ij + Îľkl = Îľijkl ) can be exploited to remove the addition of

both terms, now written as the product of two terms,
ij

ij

ab kl cd
ab kl cd
ab
cd
1 X Vab Vij Vcd Vkl
1 X ij ab kl cd Îľij + Îľkl
1 X Vab Vij Vcd Vkl
+
=
Vab Vij Vcd Vkl  2
ab Îľabcd Îľcd
ab Îľabcd Îľab
16
16
16
Îľ
Îľ
cd
ij
ij
ij
ijkl
kl
ijkl
Îľab
Îľabcd
abcd
abcd
abcd
ij
ijkl Îľkl
ijkl
ijkl
ijkl
ij ab
cd
1 X Vab Vij 1 X Vkl
cd Vkl = hΨ(1) |Ψ(1) iE (2) .
=
¡
 2
n
n
n
4
4
Îľcd
kl
abij Îľab
cdkl
ij

(2.65)

In diagrammatic form, these sums are represented by internal lines between vertices. The
common resolvent in each term, drawn as a line through the relevant state, is removed, and

31

the common diagrams which result are shown as a product (see [18] for more details),

R0

R0

c

k

d

l

R0
a

i

b

+

j

R0

R0

R0
a

i

j
c

b

k

d

l

=

a

R0

i

b

j ×c

k

d

l

R0

. (2.66)

R0

A similar factorization can be performed on the wave function terms. The following
example uses two similar third-order terms with different time-ordered interaction vertices.
Once again, the additive property of the energy denominators is used to factor the common
denominator between the terms, resulting in the product of two lower-order terms,
ab kl cd
ab kl cd
ab
cd
1 X Vij Vcd Vkl ab
1 X ab kl cd Îľij + Îľkl
1 X Vij Vcd Vkl ab
V
V
|ÎŚab
|ÎŚ
i
+
|ÎŚ
i
=
V
2
ij cd kl 
ij i
ij
ij
cd
abcd
ab
ab
abcd
ab
16
16
16
Îľ
Îľ
Îľ
Îľ
Îľ
Îľ
abcd
cd
ab
Îľij Îľijkl Îľkl
abcd kl ijkl ij
abcd ij ijkl ij
abcd
ijkl
ijkl
ijkl
ab
(1)
Vcd
1 X Vij
1 X Vkl
|Ψn i (2)
ab
cd
kl
=
=
En
 2 |Όij i ¡
4
4
Îľn
Îľcd
kl
abij Îľab
cdkl
ij
R0
R0
R0

a

i

b

R0

j
c

k

d

l

+

R0

R0
a

i

b

j

c

R0

k

d

l

a

=

i

b

j× c

k

d

l

R0

R0

(2.67)
The factorization theorem is also valid with off-diagonal Fock terms and applies to the
MBPT expansions of both the wave function and energy. Therefore, the MBPT expansions
in Eqs. (2.61) and (2.62) can be written in terms of linked diagrams only [7],
∞ h
in
X
|Ψi =
R̂0 (V̂ − ∆E0 ) |Φ0 iL ,

(2.68)

n=0

∆E0 =

∞
X

h
in
hΦ0 | V̂ R̂0 (V̂ − ∆E0 ) |Φ0 iL ,

n=0

32

(2.69)

where “L” denotes that no diagrams with closed, disconnected pieces should be included.
This result not only simplifies the MBPT expressions, but it guarantees the size-extensivity of
the MBPT wave function at each order [18]. Also, it is a useful step towards coupled-cluster
theory which reorganizes the connected diagrams from MBPT such that certain classes can
be summed to infinite order, see section 3.2.

33

Chapter 3
Coupled-Cluster Theory
Coupled-cluster theory is a powerful method for approximating solutions to the many-body
Schrödinger equation. Because of its effectiveness and economical scaling it has been a staple
of many-body quantum mechanics for decades. This chapter details various aspects of the
coupled-cluster approach and presents ground-state results for multiple systems. First, the
CC wave operator and the corresponding effective Hamiltonian will be introduced and used
to derive the CC equations. Then, results from MBPT are used to illuminate the underlying
many-body physics of the CC wave function. After the mathematical foundations of coupled
cluster theory are outlined, specific implementation details are discussed and demonstrated
by focusing on two simple examples. Finally, the nuclear many-body problem is formally
introduced, along with a brief description of the nuclear interactions used in this work, and
selected results are shown.

3.1

Exponential Ansatz

Coupled cluster theory is based on expressing the A-particle correlated wave function |Ψi
using the exponential ansatz [9, 19, 15, 16],

|Ψi = eT̂ |Φ0 i.

34

(3.1)

The cluster operator T̂ ≡ T̂ 1 + T̂ 2 + · · · + T̂ A , is composed of k-particle k-hole excitation
operators, T̂ k , which have the form,

T̂ 1 ≡
T̂ 2 ≡

X

tai

n

†
âa âi

o

,

ai

o
1 X ab n † †
tij âa âb âj âi ,
4
abij

..
.

T̂ k ≡


o
1 2 X a1 ···ak n †
†
âa1 · · · âak âi · · · âi ,
t
1
k
k! a ...a i1 ¡¡¡ik

(3.2)

1 k
i1 ...ik

a ...a

where the unknown matrix elements, ti 1...i k , are known as cluster amplitudes [18].
1 k
Using the CC ansatz, the Schrödinger equation can be rewritten by left-multiplying with
hΦ0 | e−T̂ as,
Ĥ eT̂ |Φ0 i = E eT̂ |Φ0 i
−→ hΦ0 | H |Φ0 i = E,

(3.3)

where the coupled cluster effective Hamiltonian is defined as,

H ≡ e−T̂ Ĥ eT̂ ,

(3.4)

in which the wave operator, eT̂ , acts as a similarity transform on the Hamiltonian. Using
the normal-ordered Hamiltonian HN and the correlation energy ∆E from Eq. (2.4), this can

35

rewritten as,

ĤN eT̂ |Φ0 i = ∆E eT̂ |Φ0 i
−→ hΦ0 | HN |Φ0 i = ∆E,

(3.5)

where the normal-ordered effective Hamiltonian is constructed equivalently to Eq. (3.6),

HN ≡ e−T̂ ĤN eT̂ .

(3.6)

An important characteristic of the effective Hamiltonian, H and HN , is that because
the cluster operator, which contains no de-excitations, is not Hermitian, the exponential
wave operator cannot be unitary, and thus H is not Hermitian. This initially seems like
an explicit contradiction to any standard quantum mechanics formulation where observables
are associated with the real eigenvalues of Hermitian matrices. Technically, the existence
and reality of the CC solution is only guaranteed when the full cluster operator is used
T̂ = T̂ 1 + · · · + T̂ A , see [106, 107, 108, 109]. However, while the non-Hermiticity does
require some special considerations which will be discussed, this fundamental problem does
not materialize in this work.
The effective Hamiltonian in Eq. (3.6) can be rewritten with commutators ([Â, B̂] =
ÂB̂ − B̂ Â) according to the Baker–Campbell–Hausdorff expansion as,
H = Ĥ + [Ĥ, T̂ ] +

1
1
1
[[Ĥ, T̂ ], T̂ ] + [[[Ĥ, T̂ ], T̂ ], T̂ ] + [[[[Ĥ, T̂ ], T̂ ], T̂ ], T̂ ] + · · · ,
2!
3!
4!

(3.7)

which terminates at four-nested commutators when using a two-body interaction. This
commutator expression ensures that CC theory is size-extensive and contains only connected

36

terms. In addition, because T̂ is an excitation operator, terms of the form T̂ Ĥ are necessarily
disconnected and thus vanish [18]. Therefore the CC effective Hamiltonian can be further
reduced to



H = ĤeT̂ ,
c

(3.8)

where the subscript “c” indicates that only connected terms are used.

3.1.1

The Coupled Cluster Equations

In practice, the cluster operator T̂ must be truncated for calculations to be computationally
feasible. In this work, we use only single and double excitations where applicable,

T̂ = T̂ 1 + T̂ 2 .

This is known as coupled cluster with singles and doubles (CCSD), with an asymptotic


computational cost that scales like O n4p n2h , where nh is the number of hole states and
np is the number of particle states. This truncation has been successfully applied to many
problems in quantum chemistry [110] and nuclear physics [49, 111].
The unknown cluster amplitudes in CCSD, tai and tab
ij , can be calculated by left-multiplying
−T̂ , respectively,
Eq. (3.1) with hΦai | e−T̂ and with hΦab
ij | e

hÎŚai | H |ÎŚ0 i = 0,
hÎŚab
ij | H |ÎŚ0 i = 0.

(3.9)

After the Fock matrix has been diagonalized, the diagonal components of Eq. (3.9) can be

37

separated and, after expanding the exponent in Eq. (3.8), the non-vanishing terms of the
CCSD amplitude equations in the HF basis become,




1
1
Ĥ2 T̂ 1 + T̂ 2 + T̂ 1 T̂ 2 + T̂ 21 + T̂ 31
2!
3!



(3.10)
|ÎŚ0 i = Îľai tai
c

 
1 2
1 2 1 3 1 2
1 4
ab
ab
hΦij | Ĥ2 1 + T̂1 + T̂2 + T̂1 + T̂1 T̂2 + T̂2 + T̂1 + T̂1 T̂2 + T̂1
|ÎŚ0 i = Îľab
ij tij ,
2
2!
3!
2!
4!
c
hÎŚai |

where Îľ are equivalent to the MBPT energy denominators from Eq. (2.7).
These non-linear equations are solved using an iterative procedure where the cluster
amplitudes on the right-hand side of Eq. (3.9) and Eq. (3.10) are updated by calculating
the terms on the left-hand side until a fixed point is reached. Like the HF iterative procedure, employing convergence acceleration techniques can reduce the number of CC iterations
required. Detailed techniques for solving these equations are discussed in section 3.4.

3.2

Linked-Cluster Theorem and MBPT

The exponential ansatz in Eq. (3.1) and the cluster amplitudes in Eq. (3.1) are not just useful
mathematical constructs for solving the many-body problem. They represent physical manybody dynamics and can be derived from the results of many-body perturbation theory, see
section 2.7. The linked-cluster theorem [16, 104, 105] states that the sum of all time orderings
of a disconnected diagram is equal to the product of the two subdiagrams. Using the results
and techniques from the factorization theorem, 2.7.1, the linked-cluster theorem can be used
to factorize specific MBPT terms such that they can be analytically summed to infinite
order. This infinite summation is the main principle behind coupled cluster theory and can
be shown to lead naturally to the exponential ansatz.
The first step in deriving the exponential ansatz is to group all the linked MBPT diagrams
38

of Eq. (2.68) into classes according to their number of disconnected pieces. The first class,
where all the terms are connected, can be defined as the cluster excitation operator T̂ ,
depicted by the vertex type

. The connected terms can then be characterized by their

excitation type such that T̂ k corresponds to kp -kh excitations.
The T̂ 1 operator represents the class of all connected 1p -1h MBPT diagrams, which are
determined by the number of particle-hole pairs (or pair of up- and down-lines) at the top
of the diagram,

T̂ 1 |Φ0 i ≡



(1)
(2)
= T̂ 1 + T̂ 1 + · · · |Φ0 i =

+

+

+

+

+

+ ¡¡¡ ,

+

(3.11)

where only first-order T̂ (1) and second-order T̂ (2) terms are shown while resolvant lines and
labels are removed for clarity. The T̂ 2 operator similarly represents the class of all connected
2p -2h MBPT diagrams. The first- and second-order terms that belong to this class are,

T̂ 2 |Φ0 i ≡



(1)
(2)
= T̂ 2 + T̂ 2 + · · · |Φ0 i =

+

+

+

+

+

+

+

+ ¡¡¡ .
(3.12)

Additionally, the T̂ 3 operator represents the class of all connected 3p -3h MBPT diagrams,

39

of which there is only two second-order terms,

T̂ 3 |Φ0 i ≡

=



(2)
T̂ 3

+ ¡¡¡



|ÎŚ0 i =

+ ¡¡¡ .

+

(3.13)
So far, this is merely a redefinition of the connected class of MBPT diagrams up to all orders
and is not particularly useful. But the disconnected diagrams, neglected up to this point,
can be recombined using the factorization theorem 2.7.1 to provide a powerful simplification
(see [18] for a more thorough derivation).
As an example, the remaining disconnected first- and second-order diagrams can be written as the product of connected diagrams. First, the second-order disconnected diagram from
the term fˆ2 |Φ0 i can be rewritten by adding a duplicate with the left and right subdiagrams
exchanged, which doesn’t change the value because the state labels are generic. Then, the
two disconnected diagrams can be rewritten using the factorization theorem following the
example of Eq. (2.7.1),

R0
R0

b
a

i

j

1
= 
2


R0
R0

b
a

i

j

R0

+

R0

a

i
b

j



 1
= 
 2

R0

2
a

i


 .

(3.14)

The resulting product involves the first-order component of the T̂ 1 operator in Eq. (3.11),

40

(1)

T̂ 1 . Algebraically, this process is written below,
X fia fjb
X fia fjb
X fia fjb
X
Îľb + Îľai ab
ab i = 1
ab i + 1
ab i = 1
af b j
|ÎŚ
|ÎŚ
|ÎŚ
f
i j ab a b |ÎŚij i
ab Îľa ij
ab Îľa ij
ab Îľb ij
2
2
2
Îľ
Îľ
Îľ
Îľij Îľi Îľj
abij ij i
abij ij i
abij ij j
abij


!
b
1 X fia a X fj b  1  (1) 2
=
|ÎŚ i
|ÎŚj i =
.
(3.15)
T̂ 1
2
Îľai i
2
Îľbj
ai

bj

This procedure can be repeated for the single second-order disconnected term from
V̂ 2 |Φ0 i,

R0
R0

c
a

i

b

j

k

d

l



1
= 
2

R0
R0

c
a

i

b

k

d

l

j

R0

+


1

2

=

a

i

b

c

R0

R0

j
k

d

l





2
a

i

b

j


 .


(3.16)

A similar results gives the product involving the first-order component of the T̂ 2 operator in
(1)

Eq. (3.12), T̂ 2 . Again, the factorization process is written algebraically,
ab cd
ab cd
ab cd
1 X Vij Vkl abcd
1 X Vij Vkl abcd
1 X Vij Vkl abcd
|ÎŚijkl i =
|ÎŚijkl i +
|ÎŚ
i
16
32
32
Îľabcd Îľab
Îľabcd Îľab
Îľabcd Îľcd ijkl
ij
ij
abcd ijkl
ijkl

abcd ijkl
ijkl

abcd ijkl kl
ijkl

cd
ab
1 X ab cd Îľkl + Îľij
=
Vij Vkl abcd ab cd |ÎŚabcd
ijkl i
32
Îľ
Îľij Îľ
abcd
ijkl



ab
1 X Vij

1
= 
2 4

Îľab
abij ij

ijkl

kl



|ÎŚab
ij i

!

2
1 X Vcd
kl |Φcd i = 1 T̂ (1)
.
2
4
2
Îľcd kl
cdkl

(3.17)

kl

Lastly, for the disconnected terms from V̂ fˆ|Φ0 i and fˆV̂ |Φ0 i, the first duplication step

41

can be skipped and the diagrams can be factorized following the procedure in Eq. (2.7.1),

R0
R0

c
a

i

b

j

k

R0

+

a

i

b

R0

j
c

R0

k

c

k

a

=

i

b

j

R0

.

×

(3.18)

This factorization results in a mixed term between the first-order components from the T̂ 1
and T̂ 2 operators,
ab c
ab c
c
ab
1 X Vij fk abc
1 X Vij fk abc
1 X ab c Îľk + Îľij
Vij fk abc ab c |ÎŚabc
|ÎŚijk i +
|ÎŚijk i =
ijk i
abc
ab
abc
c
4
4
4
Îľ Îľij
Îľ Îľ
Îľ Îľij Îľ
abc
ijk

ijk

abc
ijk

ijk k

!

=

ijk

abcd
ijkl

X fc
k |ÎŚc i
Îľck k
ck

1
4

k



X Vab
ij

|ÎŚab
ij i
ab
Îľ
abij ij

(1) (1)

= T̂ 1 T̂ 2

(3.19)



(1) 2
Adding these factorized contributions from Eqs. (3.15), (3.17), and (3.19) gives 21 T̂ 1
+




(1) 2
(1)
(1) (1)
(1) 2
= 12 T̂ 1 + T̂ 2
T̂ 1 T̂ 2 + 21 T̂ 2
. Repeating this procedure for all diagrams with
two disconnected parts (L = 2) adds similar terms of all orders. The final contribution from
the L = 2 class of diagrams with two disconnected parts is,
∞ h
in
X
1
R̂0 (V̂ − ∆E0 ) |Φ0 iL=2 = T̂ 2 .
2

(3.20)

n=0

A similar form results from any class of diagrams with k disconnected pieces,
∞ h
in
X
1
R̂0 (V̂ − ∆E0 ) |Φ0 iL=k = T̂ k .
k!

(3.21)

n=0

Therefore, summing over all classes of diagrams gives the final result that justifies the

42

exponential ansatz,

|Ψi =

∞ X
∞ h
X
k=0 n=0

∞
in
X
1 k
T̂ |Φ0 i = eT̂ |Φ0 i.
R̂0 (V̂ − ∆E0 ) |Φ0 iL=k =
k!

(3.22)

k=0

This equation shows the true strength and elegance of coupled cluster theory. By an ingenious
reorganization and factorization of certain MBPT diagrams, the exponential ansatz captures
the contribution of these excitations to infinite order. A more comprehensive derivation of
the linked-cluster theorem can by found in [18].

3.3

Example: Pairing Model

It’s beneficial to illustrate a simplified example of coupled cluster theory. For this purpose,
we turn our attention to the simple pairing model. This system uses a model space of N
shells, or degenerate groups of single-particle states, each with two opposite spin orbitals.

E=δ(p-1)

p=N

EF
p=A/2

p=2
p=1
σ=-1/2

σ=+1/2

Figure 3.1: Schematic representation of the pairing model space. The shells are equally
spaced and doubly degenerate with one spin-up and one spin-down state.

43

With a closed-shell reference state, the Hamiltonian is restricted to interact only between
unbroken pairs, which can be written as,

(1)Ĥ =

δ

N
X
p

(p − 1)

h

†
âp+ âp+ +

†
âp− âp−

i

, and

N

(2)Ĥ =

−

gX † †
âp âp â â ,
2 pq + − q− q+

(3.23)

where δ and g are free parameters and the ¹ labels represent the spin-up and spin-down
states, respectively.
As with all other coupled cluster calculations in this work, the first step is transforming
the problem to the Hartree-Fock basis. In this case, the restriction to unbroken pairs means
that the original single-particle states do not mix with states in other shells. In fact, the
original basis is already an eigenbasis of the Fock operator, Eq. (2.5), so that the HF transformation is reduced to a redefinition of the single-particle energies to their corresponding
Hartree-Fock energies, leaving the two-body interaction unchanged,

pmp

Îľpmp = (1)Hpmp +

X
imi

(2)H pmp imi
pmp imi

pq

g
= δ (p − 1) − ,
2

pq

Vrs = (2)Hrs .

(3.24)

Because of the pairing restriction, only hole-state energies have to be redefined.
The next step in calculating the ground-state correlation energy is to solve the CCD
equations in the Hartree-Fock basis. The system of equations comes from the terms of Eq.
(3.10) that contain only the T̂ 2 operator, and are most easily derived with diagrammatic

44

techniques, see [18].

j

hÎŚab
ij |

(Ĥ eT̂ 2 )C |Φ0 i

i

+

i

+

i
d

c

c

+

4

ab cd
Vkl
cd tkl tij + P̂ (ab)

klcd

X
klcd

k

+

b

a
d

b

i
c

k

c

j

j

l

j

l

a
b

+

a
j

i

i

=

b

k

b

a

j
c

a

+

d

ab
(εi + εj − εa − εb ) tab
ij = Vij +

1X

a

i
j

k

i

−

b

c

j

b
k

l

l

i

a
b

k

+

−

j

a

a

a

b

d

i
j

+
c

k

l

b

d

X
1 X kl ab 1 X ab cd
ac
Vij tkl +
Vcd tij − P̂ (ij|ab)
Vkb
ic tkj
2
2
kl

cd

ac bd
Vkl
cd tlj tki − P̂ (ij)

1X
2

klcd

kc

ab cd
Vkl
cd tlj tki − P̂ (ab)

1X
2

db ca
Vkl
cd tij tkl .

klcd

(3.25)
The CCD equations are written in this particular form so that an initial guess for all the
amplitudes tab
ij can be used to calculate the right-hand side of Eq. (3.3) and update the
amplitudes on the left-hand side iteratively until the amplitudes do not change within a
certain tolerance.
Lastly, the CCD correlation energy can be found with the T̂ 2 term of Eq. (3.1),

∆ECCD = hΦ0 | (Ĥ eT̂ 2 )C |Φ0 i =

c

k

d

l

=

1 X kl cd
Vcd tkl .
4

(3.26)

klcd

The correlation energies for a specific case, with δ = 1.0, N = 4, and A = 4, were calculated
for a number of different interaction strengths, g. The results are shown in Fig. 3.2 along
45

with the MBPT correlation energies to third (MBPT3) and fourth (MBPT4) orders for
comparison. The nonzero MBPT expressions for second and third order are,
ij

ab
1 X Vab Vij
,
∆EMBPT2 =
4
Îľab
ij

(3.27)

ijab

ij

ij

ab cd
kl ab
1 X Vab Vcd Vij
1 X Vab Vij Vkl
∆EMBPT3 = ∆EMBPT2 +
+
.
cd
ab Îľab
8
8
Îľab
Îľ
Îľ
ij ij
ij
ijabcd

ijklab

(3.28)

kl

Generally, the fourth-order expression contains 39 additional terms. However, most of these
vanish in this case because of the form of the pairing interaction.

Exact
MBPT3
MBPT4
CCD

Correlation energy

0.0
0.1
0.2
0.3
0.4
0.5
1.0

0.5

0.0

Interaction strength, g

0.5

1.0

Figure 3.2: Correlation energy for the pairing model with exact diagonalization, CCD, and
perturbation theory to third (MBPT3) and fourth order (MBPT4) for a range of interaction
values, g.
Also shown are the exact results from the CI method. With an example this small, it’s
possible to diagonalize, and even show explicitly, the full CI Hamiltonian matrix for an exact
result. There are six possible Slater determinants with no broken pairs, one representing the
46

reference state, four representing various 2p -2h excitations, and one representing a 4p -4h
excitation. The diagonal elements of the matrix include the single-particle energies of the
constituent states, and the matrix elements between Slater determinants with no overlap
vanish in accordance with the Slater-Condon rules, see Eq. (2.6) and [89, 99]. The full
Hamiltonian matrix to be diagonalized is,


2δ − g −g/2 −g/2 −g/2 −g/2
0
 −g/2 4δ − g −g/2 −g/2
−0
−g/2 


 −g/2 −g/2 6δ − g
0
−g/2
−g/2 


H=
.
 −g/2 −g/2
0
6δ − g −g/2
−g/2 


 −g/2
0
−g/2 −g/2 8δ − g −g/2 
0
−g/2 −g/2 −g/2 −g/2 10δ − g

(3.29)

As methods to obtain the ground-state correlation energy, both CI and CC decouple, to
some degree, the reference state and excitations from it. This decoupling has the effect of
shuffling correlations into the reference state and suppressing matrix elements connected to
it. Full decoupling between the reference state and all other Slater determinants can only
be achieved with untruncated versions of these techniques, while decoupling of the strongest
correlations can be approximately obtained with appropriate truncations. However, unlike
other many-body methods, the most unique aspect of the CC similarity transformation is
that because of its non-unitarity, the resulting Hamiltonian will not be Hermitian. This
decoupling and non-Hermiticity can be seen in Fig. 3.3, which shows the effect of the CC
similarity transformation on the Hamiltonian for a pairing case with N = 6 and A = 4.
The effective Hamiltonian H shown in Fig. 3.3 can be explicitly built, and it happens to be
beneficial to do so as part of most CC calculations, both for solving the CC equations and
for use in post-CC methods. This topic is discussed in detail in the next section.

47

4p-4h

0p-0h

4p-4h

3

3

0

0

-33

-33

-66

-66

-99

-99

4p-4h

4p-4h

2p-2h

2p-2h

2p-2h

2p-2h

0p-0h

12
-12

12
-12

log(H)

log(H)

Figure 3.3: Visualization of the CCD similarity transform on the pairing Hamiltonian for
four particles and six shells. This shows the main function of CCD, which is to decouple
2p -2h excitations from the ground state, shown by the suppression of matrix elements on the
first column. In the pairing model, this also has the effect of decoupling 2p -2h excitations
from 4p -4h excitations. Also, the non-unitary nature of the transformation is obvious given
the asymmetry of the resulting Hamiltonian.

3.4

Solving the Coupled Cluster Equations

As described above, the coupled cluster equations are solved by first initializing all of the
cluster amplitudes, then updating them by computing various sums over particle and hole
combinations, CC (T̂ ). This updating procedure is performed over multiple iterations until
the amplitudes stay unchanged within a certain tolerance,

Initialize :

T̂ (0) = T̂ init ,

Iterate :

T̂ (n+1) ← CC (T̂ (n) ) .

(3.30)

Generally, the convergence of this algorithm depends on the relative magnitudes between
the inter-particle force and the single-particle energy spacing. For a relatively small Fermi
gap between the closed shell of occupied states and the unoccupied particle states, the energy
denominators will approach zero and cause a divergent or chaotic solution [112]. Physically,

48

this situation occurs when a system exhibits strong many-particle clustering, which is difficult
to capture with only the single and double excitations of CCSD. A simple way to avoid this
ill-defined behavior is to scale the energy denominators for early iterations or to employ
linear mixing to dampen the solution,

T̂ (n+1) ← αCC (T̂ (n) ) + (1 − α) T̂ (n) .

(3.31)

If a solution to a highly-collective system does not diverge, it typically converges very slowly.
To improve the convergence rate, techniques already utilized for the Hartree-Fock iterations
can also be employed here, such as DIIS [96, 97] or Broyden’s method [98]. The additional
computational complexity for the CC iterations is simply a multiplicative factor equal to the
number of iterations performed. Typical calculations in this work with DIIS acceleration
are converged within ∟ 30 iterations. Therefore, any significant improvements to the CC
algorithm will involve the expensive sums embedded within the function CC (T̂ (n) ).

3.4.1

Symmetry Channels

For the coupled cluster equations, as well as many other many-body methods, the first way
to simplify the various sums is to exploit any symmetries of the underlying Hamiltonian of a
particular problem. These symmetries manifest as conserved quantities, and because of the
underlying nature of the cluster operators, see Section 3.2, these must conserve these quantum numbers as well. For example, the pairing Hamiltonian of Section 3.3 has a symmetry
that conserves both the total spin projection and the number of pairs. The Coulomb Hamiltonian of Section 3.5, which has translational symmetry, conserves the linear momentum of
any state. Finally, the spherical symmetry of the nuclear Hamiltonian ensures that angular

49

momentum and parity are conserved. To utilize these symmetries, any sums that contain
many-body states with different conserved quantum numbers can be ignored. For efficiency,
these symmetry groups can be pre-sorted into channels, ΣΞ~, where Ξ~ represents the relevant
quantum numbers of a certain channel.
For CCSD calculations, useful types of channels include the direct two-body channels,
Σξ~ –which categorizes the vector sum of two single-particle-state quantum numbers–and the
1

cross two-body channels, Σξ~ –which categorizes the vector difference of two single-particle2
state quantum numbers or, equivalently, the vector sum of a the quantum numbers of a
single-particle state and a time-reversed single-particle state.

Ξ~pq = Ξ~p + Ξ~q

−→

ξ~pq̄ = ξ~p − ξ~q = ξ~p + ξ~q̄

|pqi ∈ ΣΞ~ =Ξ~
1 pq
−→

|pq̄i ∈ Σξ~ =ξ~
2 pq̄

(3.32)
(3.33)

Also useful are the one-body channels, ΣΞ~ , which categorize both single-particle states by
3
their conserved quantum numbers. These one-body channels can also characterize a special
type of three-body state: the vector difference between the quantum numbers of a direct
two-body state and a single-particle state or, equivalently, the vector sum of the quantum
numbers of a two-body direct state and a time-reversed single-particle state,

ξ~p −→

|pi ∈ ΣΞ~ =Ξ~
3 p

ξ~pqr̄ = ξ~p + ξ~q − ξ~r = ξ~p + ξ~q + ξ~r̄ −→

(3.34)
|pqr̄i ∈ Σξ~ =ξ~ .
3 pqr̄

(3.35)

Using these channel structures, the interaction matrix elements and cluster amplitudes
can be built in different ways. The full applicability of these structures are shown in detail

50

in appendix C, but a few examples using sums in the CCD equations (3.3) are shown here.
The direct two-body channels can be used when two summed indices appear in the bra- or
ket-state of multiple matrix-elements,

1 X ab cd 1 X ab cd
Vcd tij =
Vcd tij for |cdi ∈ ΣΞ~ = ΣΞ~ .
2
2
ij
ab

(3.36)

|cdi

cd

The cross two-body channels are used when two summed indices appear in the opposite
corresponding bra- and ket-states of multiple matrix-matrix elements,

X
kc

ac
Vkb
ic tkj =

aj̄

for |kc̄i ∈ Σξ~ = Σξ~ .
Vkc̄
t
ib̄ kc̄
aj̄
ib̄
|kc̄i

X

(3.37)

Lastly, the one- and three-body channels are used when a single summed index appears
opposite a direct two-body state in multiple matrix elements,

1 X kl db ca 1 X klc̄ d a
Vcd tij tkl =
Vd tij b̄ tklc̄ for |klc̄i, |di ∈ Σξ~ = Σξ~ .
a
2
2
ij b̄
klcd

3.4.2

(3.38)

|klc̄i
|di

Matrix Structures and Intermediates

Symmetry channels not only provide an organized structure for the interaction matrix elements and cluster amplitudes, and remove any terms that violate the underlying symmetry
of a problem, but they also naturally provide an efficient way of performing sums using
matrix-matrix multiplications. For example, the sums in Eqs. (3.36)–(3.38) can be reformulated as matrix-matrix multiplications by structuring the channel-separated interaction
matrix elements and cluster amplitudes into individual matrices. These operations can be
performed very quickly using highly optimized linear algebra algorithms like those found in

51

BLAS (Basic Linear Algebra Subprograms) [113]. The matrices can be reordered so that the
summed indices correspond to the internal columns and rows of those matrices.
For the direct two-body case of Eq. (3.36), the structures are already in the correct order
such that the state |cdi, indexed by the columns of V and the rows of t, is summed by
multiplying the two matrices,

1 X ab cd 1 ab cd
Vcd tij = Vcd ¡ tij for |abi, |iji, |cdi ∈ ΣΞ~ .
2
2
1

(3.39)

cd

For the cross two-body case of Eq. (3.37), the states |ji, |bi, and |ci are time reversed so
that the summed variables are collected in a state |kc̄i. Then the matrix structures are
reordered so this state is indexed by columns and rows of t and V, respectively,

X
kc

aj̄

ac
kc̄
Vkb
ic tkj = tkc̄ · Vib̄ for |aj̄i, |ib̄i, |kc̄i ∈ Σξ~ .
2

(3.40)

Lastly, for the case of Eq. (3.38), the states |bi and |ci are time-reversed so that the states
|di and |klc̄i appear in two different matrix elements. Then, the matrix structures are
reorganized so the summed states occur in the appropriate rows and columns for matrixmatrix multiplication,

1 X kl db ca 1 a
· tdij b̄ for |ai, |ij b̄i, |klc̄i, |di ∈ Σξ~ .
Vcd tij tkl = tklc̄ · Vklc̄
d
2
2
3

(3.41)

klcd

These sums correspond to different components of the updated cluster amplitudes according to Eq. (3.30), so that different channel structures of the matrix-matrix multiplications
correspond to different amplitude structures according to the sum’s external indices. The
two external direct two-body states of Eq. (3.39), |abi and |iji, naturally map to the direct

52

amplitude structure,

tab
ij ←

1 ab cd
V ¡ t for |abi, |iji, |cdi ∈ ΣΞ~ .
2 cd ij
1

(3.42)

Similarly, the two external cross two-body states of Eq. (3.40), |aj̄i and |ib̄i, naturally map
to the cross amplitude structure,

t

aj̄
ib̄

aj̄

← tkc̄ · Vkc̄
for |aj̄i, |ib̄i, |kc̄i ∈ Σξ~ .
ib̄
2

(3.43)

Lastly, the one- and three-body external states of of Eq. (3.41), |ai and |ib̄i, naturally map
to the one-body amplitude structure characterized by the index a,

taij b̄ ←

1 a
tklc̄ · Vdklc̄ · tdij b̄ for |ai, |ij b̄i, |klc̄i, |di ∈ Σξ~ .
2
3

(3.44)

The last summation in the matrix-matrix form of Eqs. (3.41) and (3.44) involves two
multiplications, which suggests the need for an intermediate matrix to hold the result of
the first operation. This is the last main ingredient to an efficient CC algorithm. To see
the benefit of intermediate structures, it’s helpful to examine an expensive sum from the
CCD equations. For typical calculations, particle states outnumber hole states by an order
of magnitude, np ∟ 10nh , which means that one of the most expensive sums is,
1 X kl ab cd
Vcd tkl tij .
4

(3.45)

klcd

Because this term must be computed for each tab
ij , its computational cost naively scales as

53



O Nh4 Np4 . However, using the matrix form of this sum and an intermediate matrix,
1 X kl ab cd 1 ab  kl cd  1 ab
ab
Vcd tkl tij = tkl ¡ Vcd ¡ tij = tkl ¡ Xkl
ij → tij .
4
4
4

(3.46)

klcd



this term is now computed as the combination of two sums, each scaling as O Nh4 Np2 .
These intermediates can also be used as a way to combine similar sums. For example, the
last step of Eq. (3.46) has a very similar structure to the first sum in Eq. (3.3). Therefore,
the two sums can be written with a common intermediate as,


1 X kl ab 1 X kl ab cd 1 ab
1 kl cd
1
kl
ab
Vij tkl +
Vcd tkl tij = tkl ¡ Vij + Vcd ¡ tij = tab
¡ Xkl
ij → tij ,
2
4
2
2
4 kl
kl

klcd

1 kl cd
kl
where, Xkl
ij = Vij + 2 Vcd ¡ tij .

(3.47)

kl is equivalent to the hhhh
It just so happens that this form of the intermediate Xij

component of the CCD similarity transformed Hamiltonian, H . Constructing other intermediates in this way gives similar results, so it’s a natural extension to actually construct
the effective Hamiltonian at each iteration for the express purpose of using it as different
intermediate components for the CC equations. This has the added benefit of having already
computed the effective Hamiltonian for post-CC methods. The different components of the


T̂
2
CCD effective Hamiltonian, H CCD = Ĥ e
, are written below in both algebraic and
c

diagrammatic form. One-body components correspond to the vertex type
body terms correspond to the vertex type

and two-

. The pp, one-body component of H CCD

54

is,

a

a

a

=
b

+
b

k

c

l

b

1 X kl ac
Xba = fba −
Vbc tkl .
2

(3.48)

klc

The hh, one-body component is,

j

j

j

=
i

+
i

c

d

k

i

1 X ik cd
Xji = fji +
Vcd tjk .
2

(3.49)

kcd

The hhhh, two-body component, which appears as the intermediate in Eqs. (3.46) and (3.4.2),
is,

k

l

k

=
i

j

ij

k

l

+
i

ij

Xkl = Vkl +

j

c
i

1 X ij cd
Vcd tkl .
2
cd

55

l
d
j

(3.50)

Lastly, the hphp, two-body component is,

j

a

j

j

a

=
i

b

i

a

1
+
2

b

c

k

i

b

1 X ik ca
Vcb tjk .
2

ia = Via −
Xjb
jb

(3.51)

kc

Using these terms, the CCD equations can be written in pseudo-linear form using the
2p -2h component of effective Hamiltonian form of the equations, Eq. (3.9). This also explicitly shows the decoupling of the effective Hamiltonian with 2p -2h excitations. The pphh,
two-body component, which should vanish when the CCS amplitudes have converged, is,

a

i

b

j

a

i

b

a

j

=0=

i
i

+

b

j

a

+

c

a

b

i

+

i
j

c

+

ab = 0 = Vab + P̂ (ab)
Xij
ij

c

+

1X
2

cd

cd
X 0ab
cd tij +

2

kl

+

b

j
k

X

b

a

l

Xca tcb
ij − P̂ (ij)

1X

i
b

k

X

j

a

d

j

k

c

Xik tab
kj

k

kl tab − P̂ (ab|ij)
Xij
kl

X

kb tac .
Xic
kj

(3.52)

kc

The components and intermediates of the CCSD effective Hamiltonian are much more complicated and are shown with their corresponding sums in appendix B.

56

3.5

Example: Homogeneous Electron Gas

Another relatively simple calculation using the CCD approximation is the homogeneous
electron gas. This example aims to calculate the ground state energy of a three-dimensional
gas of electrons subject to Coulomb repulsion. This is an approximate model of the valence
electrons in a metal, subject to a uniform background of positive charge from the nuclei and
core electrons [114]. As will be explained below, this calculation employs pure-momentum
eigenstates such that 1p -1h excitations from the reference state are forbidden by momentum
conservation. This means that the problem reduces to the doubles approximation. To obtain
realistic results, a sufficiently-sized basis with a sufficient number of electrons must be used.
Therefore, the improvements discussed in Section 3.4 are necessary to keep the computation
time manageable as the system size increases.
With a uniform background potential, the electron gas can be constructed using eigen2

2
functions of the kinetic energy operator, (1)Ĥ = T̂ = −~
2m ∇ . In an infinite volume, however,

there are an unlistable number of plane wave modes which satisfy this condition due to the
continuous nature of the linear momentum eigenstates. Therefore, the single-particle orbits
will be constructed in a finite box of volume Ω and length L, and then the limit L → ∞ can
be taken after various expectation values have been computed,
−~2 2
∇ φkσ (r) = k φkσ (r),
2m
1
φkσ (r) = √ exp (ikr)ξσ ,
Ω

(3.53)

where m is the electron mass, k is the wave number, and Ξσ is the spin function for either

57

spin up or down electrons

Ξσ=+1/2 =

1
0

!
ξσ=−1/2 =

0
1

!
.

(3.54)

Assuming the single-particle orbits follow periodic boundary conditions within the containing box (φ(ri ) = φ(ri + L) for i = x, y, z) the wave numbers are quantized,

ki =

2πni
L

i = x, y, z

ni = 0, Âą1, Âą2, . . .

(3.55)

A state can therefore be characterized by the quantum numbers nx , nx , and nx as well as
the spin quantum number σ. The energy of such a state, independent of the spin, can be
written as

~2  2
~2
nx ,ny ,nz =
kx + ky2 + kz2 =
2m
2m





2π 2  2
nx + n2y + n2z .
L

(3.56)

Now that the single-particle orbits are established, a particular basis consisting of these
orbits can be chosen such that all states are included up to a closed shell. This basis is then
filled with electrons until a closed Fermi level is obtained. Additionally, only the unpolarized
case, in which all orbitals are occupied with one spin-up and one spin-down electron, will be
considered here. For this spherical-type level structure, the number of electrons required for
closed shells increases quickly. For example, the first six shells contain 2, 14, 38, 54, 66 and
114 states, respectively.
A finite number of electrons A in a finite box of volume Ω naturally leads to the characterization of an infinite system by its number density density ρ = A/Ω. The average

58

y

L

ky

L
Fourier
Transform

L
x

kF
kx

2π
L

z

Figure 3.4: Visulization of the Fourier transform of a finite box. This transformation characterizes the construction of the single-particle basis for infinite matter, mapping plane waves
in coordinate space onto finitely-spaced points in momentum space.
inter-electron distance, or Wigner-Seitz radius, is defined as

4 3 1
πr = ,
3 s
ρ


rs =


3 1/3
.
4πρ

(3.57)

In practice, these calculations are defined by the total number of shells included in the basis,
the number of electrons, and the Wigner-Seitz radius, usually given in units of the Bohr
~ , where c is the speed of light and Îą is the fine-structure constant.
radius, rb = mcÎą

The last ingredient to this many-body calculation is the interaction between the electrons,
the well-known Coulomb force. Using atomic units, where the elementary charge e = 1 and
1 = 1, this potential is simply
the Coulomb constant 4πε
0

V (r1 , r2 ) =

1
.
|r1 − r2 |

(3.58)

As mentioned in chapter 2, this potential can be utilized in second-quantization by computing

59

antisymmetrized integrals over the basis states. In this case, the integrals have the form,

(2)H pq
rs

≡

Z

dr1 dr2 φ∗kp σp (r1 )φ∗kq σq (r2 )



1
φkr σr (r1 )φks σs (r2 ) − φks σs (r1 )φkr σr (r2 ) .
|r1 − r2 |
(3.59)

The symmetries of the Coulomb potential guarantee that the total linear momentum and
total spin projection are conserved such that,

kp + kq = kr + ks

and

σp + σq = σr + σs .

(3.60)

The integral is relatively simple given the form of the basis functions. The result is given in
terms of the momentum transfer, q1 = kp − kr and q2 = kp − kr ,
(2)H pq
rs



δσp σs δσq σr
4π~cα δσp σr δσq σs
−
=
Ω
|q1 |2
|q2 |2

(3.61)

The last preparation step before performing the coupled cluster algorithm is the HartreeFock transformation. As with the pairing model, the single-particle orbitals are already
eigenfunctions of the Fock operator, in this case because the translational invariance of
pi
the plane wave basis functions ensure that the HF terms of the form (2)Hqi vanish due to

momentum conservation, see Eq. (2.5). Therefore, the HF transformation consists simply of
redefining the single-particle energies while the two-body interaction is left unchanged.

Îľp = kp +

X

(2)H pi
pi

i
pq

pq

Vrs = (2)Hrs

(3.62)

As mentioned above, in the plane-wave basis, any single excitation from the reference

60

state vanishes automatically due to momentum conservation so that tai = 0. Therefore,
it’s necessary to include only double excitations (before adding triples, etc.). Therefore,
calculations for the electron gas use the pseudo-linear form of the CCD equations (B.21)
and an effective Hamiltonian that excludes single excitations in Eqns. (3.48)-(B.4). To
explore the HEG equation-of-state, the total energy per electron can be calculated as a
function of the Wigner-Seitz radius. In the limit N, L → ∞, the plot in Fig. 3.5 represents

0.20

N = 14
N = 54
N = 246
VMC

0.15

E/A (Ha)

0.10
0.05
0.00
0.05
0.10

2

4

6

8

10

Wigner − Seitz Radius (r s )

12

14

16

Figure 3.5: CCD energy per electron in Hartrees for the 3D homogeneous electron gas as
function of the Wigner-Seitz radius in units of Bohr radii. The calculation used periodic
boundary conditions and a basis with 25 shells, resulting in a total of 1238 single-particle
states. Also plotted are the variational quantum Monte Carlo (VMC) results from [3].

the equation-of-state for a 3D electron gas at absolute zero. This curve can reveal many
thermodynamic properties of the electron gas including the saturation density and saturation
61

energy, which occur at the lowest point on the curve. The CCD results are compared
with the quasi-exact results from variational quantum Monte Carlo calculations from [3].
The discrepancies between the saturation energies from the two methods can be partially
attributed to an insufficient basis size. However, even an appropriate extrapolation to an
infinite basis won’t be able to recover all of the required correlations, which suggests that
CCSDT might be necessary. Regardless of the value to the saturation energy, these CCD
results do qualitatively reproduce the saturation radius at rs ≈ 5.0.

3.6

Coupled Cluster for Finite Nuclei

The main purpose of this work is to calculate properties of atomic nuclei with coupled cluster
theory. From a many-body perspective, the main process is computing the converged cluster
amplitudes and thus the important correlations of the system. These amplitudes comprise
the CC similarity transformation, which is versatile for constructing any effective operator,
such as the Hamiltonian, that can act on the correlated system. Therefore, the first step in
calculating beta-decay properties of nuclei is solving for the ground-state wave function of
specific closed-shell nuclei.

3.6.1

Harmonic Oscillator Basis

Calculations of finite nuclei follow the basic structure of the algorithms used for the pairing
model and the homogeneous electron gas, but they also differ in some significant ways. Like
the other examples, the first step is to construct a proper single-particle basis and reference
state. Because the nuclear Hamiltonian conserves angular momentum and parity, it’s useful
to construct orbits that are eigenfunctions of these operators. For a system with no external

62

potential, like the electron gas, this suggests a basis made of plane waves. However, plane
waves do not represent the bound states of a nucleus very well. This property can be satisfied
by introducing a fictitious external potential that mimics the mean field from the collection
of nucleons. Some phenomenological potentials, like the Woods-Saxon potential, properly
consider the resonance and continuum states of realistic nuclei in addition to the bound
states. Many-body techniques discussed in this work have been applied to model spaces that
include all three types of single-particle states [115, 116] with some success. However, for
the many-body states considered in this work, it’s sufficient to consider only bound singleparticle states. Therefore, the nuclear basis will be constructed from the isotropic harmonic
oscillator,
1
V (r) = mω 2 r2 .
2

(3.63)

An eigenstate of the harmonic oscillator potential is defined by its principal quantum
number n and its orbital angular momentum quantum number l, which is denoted by the
letters s, p, d, f... for the values l = 0, 1, 2, 3... respectively. Because of spin-orbit terms in
the nuclear Hamiltonian, the orbital angular momentum is coupled to a particle spin to a
total angular momentum of j = |l + s|, which results in a degeneracy of 2j + 1 for each orbit.
This basis does not provide any simplification to eliminate single excitations, so CCSD will
be used for all the following calculations. A schematic version of this single-particle basis is
shown in Fig. 3.6. The shell structure of this basis is characterized by the energy quantum
numbers, e = 2n + l, of the HO single-particle spectrum. This can be used to define the
maximum-energy shell and the size of a HO basis with the parameter emax .
One issue with this construction, is that while the single-particle orbits are eigenstates of
the angular momentum operator and localized to the external potential, they are not translationally invariant, which is required by the nuclear Hamiltonian. This means that there
63

Protons

Neutrons
1f7/2

1d3/2
2s1/2
1d5/2

EF

1p1/2
1p3/2

EF

n

p

1s1/2

Figure 3.6: A schematic illustration of the harmonic oscillator basis used for calculations of
nuclei. Shown is an example of a initial reference state for carbon-14, with 6 protons filled to
the p3/2-subshell closure and 8 neutrons filled to the p1/2-shell closure. See text for details
on the single-particle states.
is a fictitious center-of-mass kinetic energy which must be removed from the Hamiltonian.
The COM kinetic energy can be written as the sum of one- and two-body pieces,
A
A
A
X
X
pp ¡ pq X p2p
pp ¡ pq
Pcm
T̂ cm =
=
=
+
.
2mA
2mA
2mA p<q mA
pq
p

(3.64)

The one-body piece is just a scaled form of the original kinetic energy operator, and the
two-body piece is given in a similar form to the original two-body Hamiltonian. Both can
be integrated into matrix elements like Eq. (2.2),
p2
dr1 φ∗p (r1 ) 1 φq (r1 ) ,
2mA
Z
p ¡p
(2)T pq ≡
dr1 dr2 φ∗p (r1 ) φ∗q (r2 ) 1 2 [φr (r1 ) φs (r2 ) − φs (r1 ) φr (r2 )] .
rs
mA
(1)T p
q

≡

Z

(3.65)

Subtracting the COM kinetic energy results in the intrinsic Hamiltonian for finite nuclear

64

systems,


1 X (1) p †
1 X (2) pq (2) pq  † †
Hrs − Trs âp âq âs âr
Ĥin = 1 −
Tq âp âq +
A pq
4 pqrs
1 X (3) pqr † † †
Hstu âp âq âr âu ât âs + · · · ,
+
36

(3.66)

pqrstu

This form of the bare Hamiltonian (up to the three-body force) is used in the Hartree-Fock
transformation. Then, after normal-ordering, the three-body piece is discarded, which is
referred to as a NN+3N-induced interaction. The use of a localized external potential has
further complications involving the COM wave function that are discussed in section 3.7.
A special property of this single-particle basis that can be exploited to reduce the computational complexity of the problem is the degeneracy of each orbital, due to the angular
momentum projection of each single-particle state, mj = {−j, −j + 1, · · · , j − 1, j}. According to the Wigner-Eckart theorem [117, 118], the geometrical component of a wave function,
dependent on its projection mj , can be isolated as a Clebsch-Gordon coefficient. Because
these coefficients have compact summation rules, any diagram and corresponding sum can be
written in terms of the j-orbitals instead of the single-particle mj states, commonly known as
J-scheme and M -scheme, respectively. Calculations in J-scheme require complicated angular momentum coupling, detailed in appendix D, but involve roughly an order of magnitude
fewer states compared with an M -scheme calculation in the same model space.

3.6.2

The Nuclear Interaction

Perhaps the most important component in nuclear structure calculations, and also perhaps
the most easily overlooked component from a many-body perspective, is the nuclear Hamiltonian. Further complicated by the composite nature of protons and neutrons, bound by

65

gluon exchange within the nucleon, the inter-nucleon interaction is a residual force of virtual
pion exchanges and other, more exotic processes. Early ab initio calculations avoided this
complexity by using phenomenological interactions, tuned to reproduce certain properties of
a nucleus. These phenomenological forces were effectively used for calculations using shellmodel CI and density-functional theory, but were restricted by the conditions of the fitted
parameters.
These problems, along with the success of quantum field theories in high-energy physics,
motivated the effort to describe the inter-nucleon interaction in terms of the underlying
theory of the strong force, quantum chromodynamics (QCD) [119, 120]. However, while
calculations of nuclei in terms of their constituent quarks using lattice QCD have made
some progress with increases in computing power, they have been confined to few nucleon
systems [121]. The problem is finding a way to express the high-energy QCD interactions as
low-energy forces between nucleons. Such a problem, containing two vastly different scales,
can be rewritten as an effective theory.
Chiral effective field theory (χEFT), which exploits the large difference in scales between
the low-energy regime of nuclear physics and the high-energy regime of QCD, is built from a
general Lagrangian consistent with the broken chiral symmetry of QCD [23, 4]. This broken
symmetry, a consequence of non-zero quark masses, results in several hadronic structures
including protons, neutrons, and mesons, the lightest of which is the pion, mπ ≈ 140MeV/c2
[122]. This can be exploited by systematically writing a Lagrangian as the sum of pion
exchanges of increasing order. Additional contact interactions, which represent exchanges of
heavier mesons, are also included and must be fit to low-energy nuclear data. The hierarchy of
χEFT terms, which contain 3N and higher many-body forces, are ordered by power counting
the expansion term (mπ /Λ), where Λ is an energy cutoff between the low- and high-energy

66

Figure 3.7: Diagrammatic form of the chiral EFT expansion up to N3 LO. The solid lines
represent nucleons and the dashed lines represent pions. The different vertices represent
higher-order interactions. Figure taken from [4].
scales, and is shown up to N3 LO in Fig. 3.7.
This work exclusively employs the NN force of the N3 LO interaction from Entem and
Machleidt with a cutoff of Λ = 500 MeV [123]. For most calculations, this interaction
is coupled with the N2 LO 3N interaction from Navrátil with a cutoff of Λ = 400 MeV
[124]. This NN+3N(400) interaction is successful at reproducing low- and medium-mass
nuclei, but begins to overbind beyond the sd-shell. As mentioned in the introduction, these
bare Hamiltonians exhibit strong repulsion at short ranges among high-momentum states.
Therefore to soften the interaction, the similarity renormalization group method is used to
integrate high-momentum modes out of the interaction while preserving observables [25, 26].

67

3.6.3

Ground-State Results for Nuclei

The main object of this section is to demonstrate the validity of all the ingredients which have
been discussed so far: the harmonic oscillator basis, the NN+3N(400) chiral interaction, and
the J-scheme CCSD algorithm. To accomplish this, calculations for different nuclei will be
compared to the corresponding experimental values. Additionally, results for different input
parameters will be presented to verify that the observables are independent of non-physical
variables. Once again, all results are computed with a HF-optimized basis, see section 2.5.
First, the ground-state energies should be independent of the SRG cutoff parameter
ÎťSRG . While the SRG evolution should preserve any observables, the renormalization process
induces 3N and higher-body forces which can be missed by truncations of the many-body
method in both the cluster amplitudes and the Hamiltonian. The trade-off here is that
larger cutoff parameters produce interactions that contain higher-momentum components,
which reduce a system’s convergence properties, but induce fewer many-body forces, so
that systems can be accurately described with fewer correlations. Conversely, a smaller
cutoff parameter means that solutions can be more easily converged, but it also requires a
many-body method that includes more correlations or higher-order forces [125]. To show
this effect, the ground state for oxygen-16 is shown for different cutoff parameters and for
both the NN and NN+3N-induced interactions. Accounting for both the small dependence
on the SRG cutoff parameter and the minor inaccuracies from the truncations made to the
cluster operator and the Hamiltonian, the rest of this work will use the NN+3N(400)-induced
interaction with an SRG cutoff parameter of λSRG = 2.0 fm−1 .
Next, any nuclear observables calculated with this framework should be independent of
the fictitious confining potential. This can be verified by showing the ground-state energies
of various nuclei as a function of the underlying harmonic oscillator energy, ~ω. Convergence

68

153.0

E (MeV)

153.5

-169.0
NN

λSRG = 1. 88 fm −1

-169.1

154.0

10
12

-169.3

155.0
155.5
126.0

E (MeV)

λSRG = 2. 24 fm −1

-169.2

154.5

126.5

e max

NN

-169.4
-128.7
NN + 3N
λSRG = 1. 88 fm −1

NN + 3N
λSRG = 2. 24 fm −1

-128.9

127.0
-129.1

127.5

128.0
-129.3
14 16 18 20 22 24 26
14 16 18 20 22 24 26

ω (MeV)

ω (MeV)

Figure 3.8: Ground-state energies for 16 O for the EM N3 LO NN only interaction and with
the added 3N interaction from Navrátil, both SRG softened with λSRG = 1.88, 2.24 fm−1 .
The energies are plotted for emax = 10, 12. The most obvious difference is between the NN
and NN+3N calculations, showing the importance of including 3N forces. The differences
between the cutoff parameters are resolved within ∟ 1% with the inclusion of 3N forces
and can be rectified further by including additional correlations or full 3N forces. The
experimental binding energy is shown with the grey dashed line.
is reached by increasing the size of the model space until the resulting curve is flat. Figure
3.9 shows the convergence for the doubly-magic, N = Z nuclei, 4 He, 16 O, 20 Ca, and 56 Ni,
where both protons and neutrons fill the same major shell closure. All the results converge
to a variance of < 1% at emax = 12 for intermediate values of ~ω. While a larger model

69

25.0
25.5

116
118
120
122
124
126
128
130

4 He

E (MeV)

26.0
26.5
27.0
27.5
28.0
28.5
330

440

40 Ca

8
10
12

56 Ni

460

340

E (MeV)

e max

16 O

480

350

500

360

520

370
10

15

20

ω (MeV)

25

30

540

10

15

20

ω (MeV)

25

30

Figure 3.9: Ground-state energies for doubly magic nuclei as a function of the harmonic
oscillator energy ~ω with the NN+3N(400) interaction, SRG softened with λSRG = 2fm−1 .
The energies are plotted for emax = 8, 10, 12, showing the convergence as the model space
increases. The results are independent of the underlying oscillator frequency to ∟ 1% for
emax = 12. The grey dashed line is the experimental binding energy. The overbinding of
this interaction becomes apparent as the system size increases.
space is always desirable, this level of variance justifies the use of emax = 12 for post-CC
calculations. Additionally, these results show the limitations of the NN+3N(400) interaction,
as overbinding increases with the system size, where the ground-state energies of 20 Ca and
56 Ni

differ from their experimental binding energies by ∟ 8% and ∟ 13%, respectively.

The ground-state results are also shown for singly-magic nuclei, where either the protons
70

or neutrons fill a sub-shell closure. This has the potential complication of a vanishing energy
gap between the hole and particle states, like the picture in Fig. 3.6, which causes undefined
behavior in the CC algorithm (see section 3.4). However, the subshell orbitals repel each
other when transformed during the Hartree-Fock algorithm [126], so these systems are valid
in some cases. The ground-state energies for 14 C, 22 O, and 34 Si are plotted as a function
of the underlying oscillator potential in Fig. 3.10. The smaller energy gap involved in these
systems results in stronger excitations missed by the CCSD approximation, causing further
deviations from the experimental values.
135

E (MeV)

90

14 C

140

22 O

145

95

150
100
105

155
160
10

15

20

ω (MeV)

25

30

10

15

20

ω (MeV)

25

30

260
265
270
275
280
285
290
295
300

e max

34 Si

10

8
10
12

15

20

ω (MeV)

25

30

Figure 3.10: Ground-state energies for singly magic nuclei as a function of the harmonic
oscillator energy ~ω with the NN+3N(400) interaction, SRG softened with λSRG = 2fm−1 .
The energies are plotted for different emax . The results are independent of the underlying
oscillator frequency to ∟ 1% for emax = 12. The grey dashed line is the experimental binding
energy. These results underbind with respect to their doubly-magic counterparts in Fig. 3.9.

3.7

Ground-State Center-of-Mass Factorization

While an intrinsic Hamiltonian can be built by removing the center-of-mass (COM) kinetic
energy, Eq. (3.66), there is still an inconsistency between the translational invariance of the
underlying harmonic oscillator basis and translationally-invariant nuclear many-body states
[127, 128]. This inconsistency can materialize in certain calculations in the form of spurious,

71

non-physical states.
Of course, this problem can be avoided by using more complicated basis states that obey
translational invariance, such as the use of Jacobi coordinates, but such methods are limited
to few-body problems [129, 130]. Another possible solution is to use the harmonic oscillator
basis in an untruncated space of Slater determinants up to a certain harmonic oscillator shell.
Known as the Nmax space, this treatment can be successfully applied within no-core shell
model calculations [28]. However, the factorial scaling of this method restricts its use to light
nuclei. It can be shown that in the Nmax space, the eigenstates of the intrisic Hamiltonian
are also eigenstates of the COM Hamiltonian, and any state perfectly factorizes into a COM
component and a translationally-invariant, intrinsic component,

|Ψi = |Ψin i|Ψcm i.

(3.67)

This factorization results in a compound energy spectrum, where the intrinsic component
of the spectrum is degenerate for each COM excitation. Therefore, the intrinsic spectrum can
be recovered by offsetting the COM Hamiltonian by the corresponding excitation energies,
such that the COM energies vanish, Ecm = 0. However, with truncated methods like CCSD,
this factorization is not guaranteed, and COM energies, no longer eigenenergies of the COM
Hamiltonian, can take on values Ecm 6= 0. In this case, the intrinsic spectrum is contaminated
with nonphysical, spurious states.
Because the specific form is irrelevant [131], the shifted COM Hamiltonian can be assumed
to take the form of a harmonic trap with a oscillator strength of ~e
ω , not necessarily equal
to the oscillator strength of the underlying basis ~e
ω,

Ĥcm (e
ω) =

1
3
Pcm
+ mAe
ω 2 Rcm − ~e
ω,
2mA 2
2
72

(3.68)

offset by the ground-state energy, 32 ~e
ω . In the Nmax space, the factorization in Eq. (3.67)
occurs regardless of ~e
ω , while in truncated methods like CCSD, the COM oscillator strength
is a free parameter which can be used to probe the level of COM contamination. If a
frequency exists such that Ecm (e
ω ) ≈ 0, then the wave function is approximately factorized,
and the COM wave function is in its ground state [132, 133].
The COM energy, Ecm , can be calculated by using a version of the Hellmann-Feynman
theorem by adding the COM Hamiltonian as a perturbation and computing the difference
quotient [134, 135],

Ecm (e
ω ) ≡ hΨ| Ĥcm |Ψi ≈


1 
hΨ| Ĥ + δ Ĥcm (e
ω ) |Ψi − hΨ| Ĥ − δ Ĥcm (e
ω ) |Ψi .
2δ

(3.69)

Because the operator Rcm depends only on the underlying single-particle basis regardless of
the COM oscillator frequency, it can be rewritten in terms of Ĥcm to find the relationship
between ω and ω
e,
1
ω
e2





1
3
3
ω = 2 Ĥcm (ω) − T̂ cm + ~ω .
Ĥcm (e
ω ) − T̂ cm + ~e
2
2
ω

(3.70)

Using the known value hΨ| T̂ cm |Ψi = 34 ~e
ω and the requirement that Ecm (e
ω ) = 0 gives
the following relation that relates the COM oscillator frequency to the underlying basis
frequency,
2
~e
ω = ~ω + Ecm (ω) ±
3

s


2
2
4
Ecm (ω) + ~ωEcm (ω).
3
3

(3.71)

The ground-state COM energies are plotted for 16 O and 40 Ca using the COM Hamiltonian with two different oscillator strengths in Fig. 3.11: that of the underlying basis, ~ω,
and one of the two solutions to Eq. (3.71), ~e
ω± . Of the two ~e
ω± , which are shown as the

73

solutions to Eq. (3.71) in the inset, one typically results in a large COM energy while the
other vanishes, which is plotted. Because the COM energies approximately vanish regardless of the underlying basis frequency, the COM wave function is in its ground state and
approximately factorized from the intrinsic nuclear wave function.
Unfortunately, when intrinsic states are coupled to COM excited-states, they can contaminate the spectrum of intrinsic states that are coupled to the COM ground-state. These
spurious states can be essentially removed from the ground-state spectrum with the LawsonGloeckner method [128]. When the proper COM oscillator strength is chosen such that the
COM ground-state energy vanishes, the COM Hamiltonian can be added to the intrinsic
Hamiltonian at an arbitrarily large scale, β, without changing the ground-state spectrum,

Ĥin → Ĥin + β Ĥcm .

(3.72)

When β is arbitrarily large, eigenenergies of intrinsic states coupled to COM excited states
will increase by the COM energy quanta, β~e
ω , such that they are removed from the range
of low-lying states of interest. The method will be used to remove spurious states from the
spectra of open-shell states in section 4.

74

60
50
40
30
20
10
0

f

ω (MeV)

16 O

5

10

15

20
ω (MeV)

25

30
ω
ω+
ω−
f
f

40 Ca

60
50
40
30
20
10
0

f

ω (MeV)

E CoM (MeV)
E CoM (MeV)

4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0.5

5

10

15

20
ω (MeV)

25

30
ω
ω+
ω−
f
f

10

15

20

ω (MeV)

25

30

Figure 3.11: Ground-state COM energies, Eq. (3.68), for 16 O and 40 Ca at various harmonic
oscillator frequencies with the NN+3N(400)-induced interaction with λSRG = 2.0 fm−1 at
emax = 12. Using the proper COM oscillator frequencies shows the approximate factorization
of Eq. (3.67).

75

Chapter 4
Equation-of-Motion Method
Once the ground state energy and the effective Hamiltonian have been calculated, any further
properties of the ground state can be calculated using the correlated wave function written
as an expansion of Slater determinants in the form of Eqn. (3.1). However, many of the
most interesting processes in nuclear physics involve excited-state properties. Additionally,
because the coupled cluster method requires a doubly closed-shell reference state, most topics
in nuclear physics that can benefit from an ab initio description are unreachable with the
standard coupled cluster approach alone.
These restrictions motivate a class of techniques known as the equation-of-motion (EOM)
methods [136]. Applied with the CC effective Hamiltonian, the equation-of-motion coupled
cluster (EOM-CC) method begins with a standard CC calculation of a closed-shell ground
state. Then, EOM target states are built onto the correlated ground-state wave function in
the same way that CI states were built from the reference state, Eqn. (2.48). Unfortunately,
because the CCSD similarity transformation only decouples the ground state from 1p -1h and
2p -2h excitations, these target states are still coupled to each other. Therefore, capturing
the relevant correlations to describe EOM states involves a CI-like diagonalization of the
effective Hamiltonian.

76

4.1

Equation-of-Motion States

The CC equations-of-motion states are built by applying a particular type of excitation
operator to the correlated ground-state wave function in the same way that the correlated
ground-state wave function was built by applying the exponential cluster operator to the
reference state,
|Ψµ i = R̂µ |Ψi = R̂µ eT̂ |Φ0 i.

(4.1)

The EOM excitation operator can be written as a linear combination of strings composed
of particle creation operators and hole annihilation operators in increasing order, R̂µ =
(0)R̂
Âľ

+ (1)R̂µ + (2)R̂µ + · · · . The form of these strings is determined by the structure of

the desired target state. For example, excited states maintain the particle number of the
reference state and ground state so that the constituent operator strings are kp -kh operators.
The A-particle excited-state EOM operator R̂A
Âľ is,

Âľ
R̂A
Âľ = r0 +

X
ai

Âľ ra
i

n
o 1X
n
o
†
µ r ab ↠↠â â + · · · ,
âa âi +
a
ij
b j i
4

(4.2)

abij

where (0)R̂µ = µ r0 represents the ground-state component of an excited state with the same
ab , ¡ ¡ ¡ , are the normal-ordered
conserved quantum numbers and the matrix elements, Âľ ria , Âľ rij

components of R̂µ .
In addition to excited states, particle-attached (PA) states can be reached by applying
strings of the form p, pph, ppphh, etc. to the ground state. These operator strings increase
the number of particles by one and gives the particle-attached equation-of-motion PA-EOM

77

operator R̂A+1
= (1)R̂A+1
+ (2)R̂A+1
+ ¡¡¡,
Âľ
Âľ
Âľ
R̂A+1
Âľ

=

X

Âľ ra

a

n
o
n o 1X
† †
†
Âľ
ab
ri âa âb âi + · · · .
âa +
2

(4.3)

abi

Lastly, particle-removed (PR) states can be reached by applying strings of the form h, hhp, hhhpp,
etc. to the ground state. These operator strings decrease the number of particles by one
and give the particle-removed equation-of-motion PR-EOM operator R̂A−1
= (1)R̂A−1
+
Âľ
Âľ
(2)R̂A−1
Âľ

+ ¡¡¡,
R̂A−1
=
Âľ

X
i

Âľr

i {âi } +

o
1 Xµ a n †
rij âa âj âi + · · · .
2

(4.4)

aij

With these target states, the Schrödinger equation can be written by applying the Hamiltonian to Eq. (4.1),

Ĥ|Ψµ i = Eµ |Ψµ i,
ĤR̂µ eT̂ |Φ0 i = Eµ R̂µ eT̂ |Φ0 i.

(4.5)

This equation can be written in terms of the effective Hamiltonian by multiplying with the
operator e−T̂ ,
e−T̂ ĤR̂µ eT̂ |Φ0 i = Eµ e−T̂ R̂µ eT̂ |Φ0 i.

(4.6)

Because both R̂µ and T̂ are excitation operators, containing only particle creation operators
and hole annihilation operators, no nonzero contractions can occur between them (see section
2.4). This means that the order of the two operators is inconsequential such that they
commute with each other (R̂µ T̂ = T̂ R̂µ ), which is also true for the exponential cluster
operator (R̂µ eT̂ = eT̂ R̂µ ). This property can be used to rewrite Eq. (4.6) in terms of the

78

effective Hamiltonian,

e−T̂ Ĥ eT̂ R̂µ |Φ0 i = Eµ e−T̂ eT̂ R̂µ |Φ0 i
H R̂µ |Φ0 i = Eµ R̂µ |Φ0 i.

(4.7)

Using the normal-ordered Hamiltonian, the equation can be rewritten as,

HN R̂µ |Φ0 i = ∆Eµ R̂µ |Φ0 i.

(4.8)

Up to this point, extending the CC method from the ground state to excited and openshell states amounts to an energy eigenvalue problem involving the normal-ordered effective
Hamiltonian and the EOM excitation operators R̂µ . The signature component of EOM
methods is reached by multiplying the normal-ordered ground-state Schrödinger equation,
Eq. (2.4), with R̂µ and subtracting the result from Eq. (4.8),






HN R̂µ − R̂µ HN |Φ0 i = ∆Eµ − ∆E R̂µ |Φ0 i.

(4.9)

h

i

The right-hand side of this equation can be rewritten as the commutator HN , R̂µ and the
energy difference can be defined as ωµ ≡ ∆Eµ − ∆E so that the fundamental EOM equation
is,
h

i
HN , R̂µ |Φ0 i = ωµ R̂µ |Φ0 i.

(4.10)

The name equation-of-motion refers to the resemblance of this equation with the commutatorbased Heisenberg representation of quantum mechanics. The objective behind this formulation is to reduce the dependence of EOM states on the ground state by removing common
terms between them. This reduction is accomplished by noticing that, like the commutator
79

h
i
between HN and R̂µ , the commutator between uncontracted terms of HN , R̂µ vanishes.
This simplifies the commutator in Eq. (4.11) to only connected terms just as the effective
Hamiltonian was simplified in Eq. (3.8),



HN R̂µ


c

|Φ0 i = ωµ R̂µ |Φ0 i.

(4.11)

Equation (4.6) constitutes a generalized eigenvalue problem which solves for the components of the EOM operator, µ r, and the energy difference, ωµ . If only the energy of excited
or open shell states are required, solving this equation is sufficient for such a task. However,
computing properties of EOM states requires the Hermitian conjugate of the EOM operator
†

R̂µ , and as encountered before, the non-Hermiticity of the effective Hamiltonian complicates
this effort.

4.2

Dual Solutions
†

Because the CC effective Hamiltonian is non-Hermitian ( HN 6= HN ), the eigenvalue problem
in Eq. (4.8) has a corresponding left-eigenvalue problem,

hΦ0 | L̂µ HN = hΦ0 | L̂µ ∆Eµ .

(4.12)

The operators L̂µ are de-excitation operators analogous to the corresponding right operators
R̂µ . The left EOM operator for particle-attached states has the form,

L̂A+1
Âľ

=

X
a

µ l {â } +
a a

o
1 Xµ i n †
lab âi âb â + · · · ,
2
abi

80

(4.13)

and the left EOM operator for particle-removed states has the form,

L̂A−1
Âľ

=

X
i

Âľ li

n
n o 1X
o
ij
† †
†
Âľ
la âi âj âa + · · · .
âi +
2

(4.14)

aij

Again, the matrix elements µ l are the normal-ordered components of L̂µ and map to a corresponding µ r. Because there is not a corresponding left eigenvalue problem for the ground
state in the form of Eq. (2.4), the EOM eigenvalue problem for L̂µ cannot be reduced to the
form of Eq. (4.11). This amounts to calculating additional non-contracted terms between
the Hamiltonian and the left EOM state and computing the energy difference directly from
ωµ ≡ ∆Eµ − ∆E.
The corresponding right and left EOM operators form a bi-orthogonal basis such that,

hΦ0 | L̂µ R̂ν |Φ0 i = δµν .

(4.15)

When this condition is fulfilled, any scaling of the right and left solutions by the factors
Îą and 1/Îą, respectively, also fulfills the condition. Therefore, because the normalization
of each solution is not determined uniquely, both solutions must be used when computing
properties of EOM states. This is accomplished by applying the bi-orthogonal solutions as
the identity operator,

1̂ =

X
Âľ

R̂µ |Φ0 ihΦ0 | L̂µ .

81

(4.16)

4.2.1

Induced Three-Body Interaction

Before solving the EOM equations, it’s beneficial to introduce relevant three-body effective
Hamiltonian terms which will be used in Eqs. (4.11) and (4.12). While the original threebody Hamiltonian was only used for the normal-ordered zero- and one-body pieces, the CC
similarity transformation induces higher-body interactions from contractions between the
Hamiltonian and the cluster operators, see Eq. (3.8). In the CCSD approximation, the fourbody interaction is the highest-order term, generated from the contraction of the Hamiltonian
with the term 12 T̂ 21 . Fortunately, the PA-EOM-CCSD and PR-EOM-CCSD methods only
include certain three-body interactions, which are shown below. The effective hpphpp threeo
n
† †
body interaction can couple two particle-attached operators of the form µ riab âa âb âi and
is generated from a term of the form V̂ T̂ 2 ,

a

b

i

a

=
c

j

c

d

iab = −
Xjcd

j

d

b

i

k

X jk
Vcd tab
ki .

(4.17)

k

Likewise, the effective hhphhp three-body interaction couples two particle-removed operators
n
o
†
Âľ
a
of the form rij âa âj âi ,

i

j

a

i

=
k

b

k

l

kla =
Xijb

b

l

82

a

c

X jk
Vcd tab
ki .
k

j

(4.18)

Because the pphpph structure in Eq. (4.17) scales as O(n2h n4p ), it can quickly overtake
a calculation’s memory allocation. Therefore, these structures are never actually built. Instead, the different components are summed as special intermediates when they are needed
during PA-EOM-CC or PR-EOM-CC calculations.

4.3

Solving the EOM equations

At this point, after solving for the CC ground-state wave function with a truncated cluster
operator, a full accounting of the remaining many-body correlations would scale factorially
like the full CI method, see section 2.6 and Fig. 2.3. Therefore it’s necessary to truncate
the EOM operators with the assumptions that the lower-order excitations will dominate
the EOM states. In this work, the particle-attached and particle-removed operators are
truncated at the 2p -1h and 1p -2h levels, respectively, which is referred to as the EOM(2)
truncation. Like other effective ab initio methods, this approximation can be systematically
improved by including higher-order terms. For the PA-EOM(2) method, the EOM operators
have the form,

R̂A+1
=
Âľ

X

Âľ ra

a

L̂A+1
=
Âľ

n o 1X
n
o
†
µ r ab ↠↠â
âa +
and
a b i
i
2
abi

X
a

µ l {â } +
a a

83

o
1 Xµ i n †
lab âi âb â .
2
abi

(4.19)

Likewise, the PR-EOM(2) method, the EOM operators have the form,

R̂A−1
Âľ

=

X
i

L̂A−1
=
Âľ

µ r {â } +
i i

X
i

o
1 Xµ a n †
and
rij âa âj âi
2
aij

n o
X ij n † † o
µ l i ↠+ 1
Âľl
a âi âj âa .
i
2

(4.20)

aij

The EOM matrix eigenvalue equation can be solved in a computationally practical way
with power-iteration methods. Traditionally, the Lanczos algorithm is used to produce the
lowest energy eigenvalues and their corresponding eigenvectors from a Hermitian matrix
[137]. In this case, with a non-Hermitian matrix, the generalized Arnoldi algorithm [138] is
used instead. These methods remove the need to build the entire matrix, hΦ0 | L̂µ HN R̂µ |Φ0 i,
to be diagonalized, instead relying on matrix-vector products performed as a step in an iterative process. In this work the iterative procefure was implemented with the numerical
software library ARPACK [139]. The matrix is simply the effective normal-ordered Hamiltonian HN and the vectors are EOM operators. Like the coupled cluster equations, this
matrix-vector product is best computed with diagrammatic techniques and are shown below.
The matrix-vector product for the right eigenvalue problem of the PA-EOM(2) method
consists of two components which can be seen clearly by left-multiplying Eq. (4.11) with the
particle-attached bra states, hÎŚa | and hÎŚab
i |, respectively. This has the effect of projecting
out the corresponding components Âľ r,



hΦa | HN R̂µ |Φ0 i = ωµ µ ra ,
c


hÎŚab
|Φ0 i = ωµ µ riab .
i | HN R̂µ
c

84

(4.21)
(4.22)

Like the coupled-cluster equations, the EOM equations are best derived with diagrammatic
techniques. For these purposes, the right EOM excitation operators are depicted by the
vertex type

. Equation (4.21) generates the following expressions and corresponding

diagrams, suppressing the state identifier Âľ for clarity,

ωra =

X

Xca rc +

X

c

Xck rkac +

1 X ak cd
Xcd rk
2
kcd

kc

a
a

a

=

a

+

+

c

c

k

.
c

d

(4.23)

k

The corresponding expressions and diagrams for Eq. (4.22) are,

ωriab =

X

ab r c + P̂ (ab)
Xci

X
c

c

a

b

i

− P̂ (ab)

X

a

b

kc

Xcb riac −

k

1 X ab cd
Xcd ri
2
cd

klcd

i

b

+

a

i

i

a

+

c

b

c

a

b
i

+
k

c

d

i
b

+
c

Xik rkab +

X
ak r cb − 1
ab cd
Vkl
Xci
cd tki rl
k
2

=

a

X

k

+
d

l

c

a
k

b

i

.

(4.24)

The vector-matrix product for the left eigenvalue problem of the PA-EOM(2) method
also consists of two components achieved in a similar fashion by right-multiplying Eq. (4.12)
with the particle-attached ket states, hÎŚa | and hÎŚab
i |, respectively, which has the effect of

85

projecting out the corresponding components Âľ l,

hΦ0 | L̂µ HN |Φa i = ∆Eµ µ la ,

(4.25)

Âľ i
hΦ0 | L̂µ HN |Φab
i i = ∆Eµ lab .

(4.26)

Because Eqns. (4.25) and (4.26) do not use commutators between HN and L̂µ , the diagrams
which describe these equations include disconnected diagrams and corresponding expressions
with matrix elements that share no indices. For these diagrams, the left EOM de-excitation
operators are depicted by the vertex type

. Equation (4.25) generates the following

expressions and diagrams, again suppressing the state identifier Âľ,

Ela =

X

lc Xac +

c

kcd

c
a

1 X k cd
lcd Xak
2

c

=

+
a

k

.
a

86

d

(4.27)

The corresponding expressions and diagrams for Eq. (4.26) are,

i = P̂ (ab) l X i +
Elab
a b

X
c

+

ci −
lc Xab

X

k X i + P̂ (ab)
lab
k

c

k

klcd

kc

c
b

i

=

i Xc
lac
b

X
X
1 X i cd
l Vki tcd
k X ci − 1
lcd Xab − P̂ (ab)
lcd
lcb
ab kl
ak
2
2
cd

a

X

k

+

a
b

i

c

a

d

b

c

+
b

a

c

+

b

i

a

i

i

b

k

+

i
a

+

b
a

+

d

l

c

.

k
a

i

b

(4.28)

i

The particle-removed equations are generated in exactly the same way. The matrixvector product for the right eigenvalue problem of the PR-EOM(2) method consists of two
components corresponding to the particle-removed bra states, hÎŚi | and hÎŚaij |, respectively,


hΦi | HN R̂µ |Φ0 i = ωµ µ ri ,
c


a.
hΦaij | HN R̂µ |Φ0 i = ωµ µ rij

(4.29)
(4.30)

c

Equation (4.29) generates the following expressions and diagrams,

ωri = −

X

Xik rk +

k

X
kc

c −
Xck rik

1 X kl c
Xic rkl
2
klc

i
i

=

i

+
k

i

+
k

87

c

.
k

l

c

(4.31)

The corresponding expressions and diagrams for Eq. (4.30) are,

a =−
ωk rij

X
k

− P̂ (ij)
i
i

j

a

ka r − P̂ (ij)
Xij
k

X
kc
j

X

a +
Xjk rik

c

k

1 X kl a
Xij rkl
2
kl

klcd

a

j

+

i

a

a

+

k

i

j

k

i

j
a

+
c

k

l

a
j

+
k

c +
Xca rij

X
d
ak r c − 1
rkl
Xci
Vkl
tca
ij
kj
cd
2

=

i

X

c

+
l

d

k

i
c

j

a

.

(4.32)

Finally, the vector-matrix product for the left eigenvalue problem of the PA-EOM(2)
method consists of two components which correspond with the particle-removed ket states,
hÎŚi | and hÎŚaij |, respectively,
hΦ0 | L̂µ HN |Φi i = ∆Eµ µ l i ,

(4.33)

ij

(4.34)

hΦ0 | L̂µ HN |Φaij i = ∆Eµ µ la .

Again, Eqns. (4.33) and (4.34) do not remove disconnected diagrams with commutators

88

between HN and L̂µ . Equation (4.33) generates the following expressions and diagrams,

Eli = −

X
k

k

klc

k

=

i

1 X kl ic
lc Xkl
2

lk Xki −

l

c

+
i

.

(4.35)

i

The corresponding expressions and diagrams for Eq. (4.34) are,

ij

j

Ek la = P̂ (ij) li Xa −
+

1X
2

cd

X

ij

lk Xka +

X
X ij
j
laik Xk
lc Xac − P̂ (ij)
c

k

k

X kj
X
ij
ij
ci − 1
lc Xak
lakl Xkl − P̂ (ij)
ldkl Vca tcd
kl
2
kc

klcd

k
i

j

a

=

j

a

k

j

k
a

j

+
i

l

+
i

c

+

i

k

+

j

a

i

a

a

j

c

+

j
i

i

+

l

d

k

.

c
i

a

j

(4.36)

a

These expressions must be computed for every iteration of the Arnoldi algorithm and
contain expensive sums that can benefit from the techniques of section 3.4. The major
difference in this case is that while the cluster operators and the Hamiltonian are scalar
operators that define the conserved quantum numbers of a symmetry channel, see section
3.4.1, the EOM operators are tensor operators that carry their own quantum numbers and
can change the symmetry of the reference state to that of the target state.

89

4.4

EOM-CC for 2D Quantum Dots

To examine some properties of the EOM-CC method, it’s useful to focus on the results
from the two-dimensional quantum dot. This system, also known as an “artificial atom”,
consists of A electrons confined within a 2D harmonic oscillator and subject to the interelectron Coulomb force. These nanoscale structures are easily tunable which makes them
relevant for probing many different quantum phenomena in experiments and theoretical
models [140, 141, 142], and because the quantum dot shell structure is similar to those of
atomic nuclei, they provide a useful system to quantify the impact of quantum effects at
different levels of correlation [143].

4.4.1

Quantum Dot Formalism

The circular quantum dot system is characterized by a harmonic-oscillator potential of the
form V (r) = mω 2 r2 /2, where ω is its angular frequency, r is the radial distance from
the center, and m is the electron mass. The single-particle states can be constructed as
eigenfunctions of the 2D quantum harmonic oscillator,



−~2 2 1
2
2
∇ + mω r φnm` σ (r) = nm` φnm` σ (r).
2m
2

(4.37)

These states are characterized by a radial quantum number n, a spin quantum number σ, and
because this systems conserves orbital angular momentum, a corresponding angular momentum projection number, m` . The energy of each single-particle state φnm` σ (r), independent
of the spin quantum projection, is given by,

nm` ≡ (2n + |m` | + 1)ω.

90

(4.38)

In addition to the spin projection σ, they are degenerate with respect to the shell index k,

k ≡ 2n + |m` |,

(4.39)

which labels each shell from zero with a total number of shells, K. The shell structure is

E/ω = k + 1

illustrated in Fig. 4.1.

+½

5

−½
0
1
2

3
1
−3

0

mℓ

ms
n

+3

Figure 4.1: The 42 lowest single-particle states (the first 5 shells) in the 2D harmonic oscillator basis. Each box represents a single-particle state arranged by m` , ms , and energy,
and the up/down arrows indicate the spin of the states. Within each column, the principal
quantum number n increases as one traverses upward.

The specific form of these basis states can be solved with the cylindrically symmetric
Fock–Darwin states Fnm` , which conserve the orbital angular momentum projection L̂z ≡
∂ . Apart from their spin component, the states are defined as [144],
−i ∂ϕ

√
1
ωRn|m | ( ωρ) × √ eim` ϕ , where
`
2π
s
2 × n! −%2 /2 m (m) 2
Rnm (%) ≡
e
% Ln (% )
(n + m)!

Fnm` (ρ, ϕ) ≡

(Îą)

and Ln

√

(4.40)

denotes the generalized Laguerre polynomial [145] of degree n and parameter Îą.

91

Like the homogeneous electron gas, (see section 3.5), the electrons in the quantum dot
interact with each other through the standard Coulomb interaction, V (r1 , r2 ) = 1/|r1 − r2 |,
which is expressed in atomic units where m = c = e = (4π0 )−1 = 1. With the form of the
single-particle wave functions, the second-quantized form of the two-body interaction can be
analytically calculated,

(2)H pq
rs

≡

Z

dr1 dr2 φ∗p (r1 )φ∗q (r2 )

1
[φr (r1 )φs (r2 ) − φs (r1 )φr (r2 )] ,
|r1 − r2 |

(4.41)

where φp (r) is shorthand for φnp m` σp (r).
p
With these ingredients in hand, the ground state energy, as well as particle-attached and
particle-removed energies, can be calculated for different values of the oscillator frequency,
ω, and the number of shells, K.

4.4.2

Quantum Dot Results

Here, CC results are shown along with results from the in-medium similarity renormalization
group (IM-SRG) method, and quasi-degenerate perturbation theory (QDPT) [146]. As for
all calculations in this work, each method applies a Hartree-Fock transformation of the basis
states from Eq. (4.37) before employing the primary many-body method. Before analyzing
the EOM results, the ground state energies calculated using CCSD and IM-SRG(2) are
shown in Fig. 4.2. Also included are results from Møller–Plesset perturbation theory to
second order (MP2), DMC [147], and full CI [148] for comparison where available.
With respect to the number of shells, both IM-SRG(2) and CCSD appear to converge
slightly faster than second order perturbation theory (MP2), mainly due to the presence
of higher order corrections in IM-SRG(2) and CCSD. However, the main conclusion that

92

=
= .

7.8

7.6

=
= .

128
126

MP2
IM-SRG(2)
CCSD
FCI
DMC

124
5

10

10.0

=
= .

20.4

12.5

311

15.0

=
= .

310

880

309

20.2
5

10

=
= .

882

10.0
12.5
15.0
(number of shells)

16

18

20

Figure 4.2: Ground state energy (in Hartrees) of quantum dots with N particles and an
oscillator frequency of ω calculated with several different methods.
can be reached from Fig. 4.2 is the overall convergence between vastly different methods, a
hallmark of properly-treated ab initio methods.
Additionally, the results for addition and removal energy calculations are summarized
in Fig. 4.3 and Fig. 4.4 respectively. Not only do these results reaffirm the validity of the
CCSD method amongst other ab initio methods, they also provide benchmark calculations
for the following CC-EOM calculations of open-shell finite nuclei.

4.5

Quality of EOM Solutions

Applying the CC similarity transformation to the Hamiltonian implicitly re-sums contributions from higher order excitations beyond the EOM(2)-CC truncation (3p -2h, 4p -3h, . . . for
PA-EOM and 2p -3h, 3p -4h, . . . for PR-EOM) into the lower-order excitations which comprise the EOM(2)-CC operators (1p -0h, 2p -1h, . . . for PA-EOM and 0p -1h, 1p -2h, . . . for
PR-EOM). This means that EOM states are weighted towards few-body excitations when

93

=
= .

2.7
2.6

=
= .

7.75

IM-SRG(2) + QDPT3
IM-SRG(2) + EOM2
CCSD + EOM2
DMC

7.50

2.5

7.25
10

10.0

12.5

15.0

(+)

5

=
= .

6.7
6.6

=
= .

18.2
18.0

6.5
5

10

=
= .

27.2
27.0

10.0
12.5
15.0
(number of shells)

16

18

20

Figure 4.3: Particle-attached energy (in Hartrees) of quantum dots with N particles and an
oscillator frequency of ω calculated with several different methods.

2.2

=
= .

2.1

7.50

=
= .

7.25

IM-SRG(2) + QDPT3
IM-SRG(2) + EOM2
CCSD + EOM2
DMC

7.00
10

10.0

12.5

15.0

( )

5

=
= .

5.30
5.25

=
= .

17.2

=
= .

26.2

17.0

5.20

26.0
5

10

10.0
12.5
15.0
(number of shells)

16

18

20

Figure 4.4: Particle-removed energy (in Hartrees) of quantum dots with N particles and an
oscillator frequency of ω calculated with several different methods.

94

compared to the corresponding state calculated with the CI method and the bare Hamiltonian.
Despite this improvement on the corresponding CI method, the EOM(2) truncation of
Eqs. (4.3) and (4.3) is not guaranteed to capture the primary components of certain collective
states. Therefore, the quality of the solution can be judged by measuring the amount of
overlap between the EOM solution and 1p -0h states for PA-EOM or 0p -1h states for PREOM. This overlap can be computed with partial norms,

n1-particle =

s
X
a

n1-hole =

sX
i

hΦa | R̂µ |Φ0 ihΦ0 | L̂µ |Φa i

=

s
X
a

|ra la |,

sX
hΦi | R̂µ |Φ0 ihΦ0 | L̂µ |Φi i =
|ri l i |.

(4.42)

(4.43)

i

If the single-particle overlap is small, the state is dominated by higher-order excitations and
requires a less-restricted truncation in order to describe it properly, which can be accomplished directly or perturbatively [149, 52]. On the other hand, large single-particle partial
norms indicate that the EOM truncation is reasonable for the relevant state.
Figure 4.5 shows this relationship using the difference between EOM-CC and full CI
energies for quantum dot particle-attached and particle-removed states, ∆E/E, as function
of the single-particle overlap. The energy difference with the exact FCI result grows as the
single-particle character of the wave function shrinks. Therefore, care should be taken when
encountering unphysical states with a small single-particle overlap.

95

0.05

PA-EOM
PR-EOM

∆E/E

0.04
0.03
0.02
0.01
0.00 0.0

0.1

0.2

0.3

0.4

Figure 4.5: Energy difference of particle-attached
EOM-CC method and the exact FCI method for
plotted against the single-particle overlap of the
(4.42) and (4.43). The strong correlation shows
judged by this metric.

4.6

0.5

n 1 − particle , n 1 − hole

0.6

0.7

0.8

and particle-removed states between the
a quantum dot with various parameters
FCI state, n1-particle or n1-hol , see Eqs.
that the quality of EOM states can be

EOM-CC for Finite Nuclei

Like its ground-state counterpart, the EOM extension to coupled cluster theory was first
applied to nuclear systems [150, 151, 152]. However the progress was again halted by the nonperturbative nucleon-nucleon interaction, so again it quickly gained prominence in quantum
chemistry [153, 18]. Since the introduction of SRG-softened chiral interactions, the EOMCC method has been established as a powerful and reliable method for reaching open-shell
systems and excited states [44, 50, 133, 154, 155]. This section shows results for the PAEOM-CC and PR-EOM-CC methods calculated with the NN+3N(400) interaction, SRG
softened with λSRG = 2.0 fm−1 . The initial closed-shell calculations were performed for 16 O
and 22 O.
The ground-state energies for particle-attached nuclei are shown in Fig. 4.6. Unlike
96

the ground states for closed-shell nuclei, which is uniquely determined, the particle-attached
ground state must be identified as the state corresponding to the lowest eigenvalue of the
converged solution. All the results are in close agreement with the experimental values while
the subshell nuclei, 23 O and 23 F, are slightly underbound like their closed-shell counterpart,
16 O

(see Fig. 3.10).

115

110

E (MeV)

17 O

120

115

125

120

130

125

135
140

130
145

E (MeV)

145

23 O

150

8
10
12

23 F

155

150

160

155

165

160

170

165
170

e max

17 F

175
10

15

20

ω (MeV)

25

30

180

10

15

20

ω (MeV)

25

30

Figure 4.6: Ground-state energies for the particle-attached nuclei 17 O, 17 F, 23 O, and 23 F
as a function of the harmonic oscillator energy ~ω with the NN+3N(400) interaction, SRG
softened with λSRG = 2fm−1 . The energies are plotted for emax = 8, 10, 12, showing the
convergence as the model space increases. The grey dashed line is the experimental binding
energy.

97

The ground-state energies for particle-removed nuclei are shown in Fig. 4.7. These
results follow the same pattern as the particle-attached states, converging to within ∟ 1%

E (MeV)

at emax = 12 and overbinding for EOM states constructed from the 22 O ground state.

100
102
104
106
108
110
112
114
116
118
110

E (MeV)

115

95
15 N

8
10
12

100
105
110
115
130
21 N

135

120

140

125

145

130

150

135

155

140

e max

15 O

10

15

20

ω (MeV)

25

30

160

21 O

10

15

20

ω (MeV)

25

30

Figure 4.7: Ground-state energies for the particle-removed nuclei 15 N, 15 O, 21 N, and 21 O as
a function of the harmonic oscillator energy ~ω with the NN+3N(400)-induced interaction,
SRG softened with λSRG = 2fm−1 . The energies are plotted for emax = 8, 10, 12, showing
the convergence as the model space increases. The grey dashed line is the experimental
binding energy.

The problem of center-of-mass contamination discussed in section 3.7 extends to the
EOM states as well. To ensure that the EOM wavefunction has been sufficiently factorized

98

from the COM wave function, an equivalent diagnostic to Fig. 3.11 can be performed. First,
the ground-state COM energy is computed using Eq. (3.69) and the corresponding COM
oscillator strength is obtained from Eq. (3.71). By using A + 1 particles for PA-EOM and
A − 1 particles for PR-EOM in the intrinsic Hamiltonian, Eq. (3.66), the A-body problem is
not properly treated. However, the approximate COM factorization for the A + 1 and A − 1
systems, according to Eq. (3.67), is reached when the COM energy vanishes at this new
oscillator strength. The results from the COM diagnostic on the PA-EOM and PR-EOM
states is shown in Fig. 4.8.
With the underlying COM oscillator strengths determined from Fig. 4.8, the excitation
spectra of the corresponding open-shell nuclei can be treated with the Lawson-Gloeckner
method [128] by artificially raising COM-coupled, spurious states out of the range of interest
by adding the COM Hamiltonian to the intrinsic Hamiltonian, Eq. (3.72). The EOM lowlying excited states of 17 O and 17 F are shown in Fig. 4.9 with and without the LawsonGloeckner term. When the term is added, the spurious negative-partiy states are removed
from the lower portion of the spectra. The non-spurious states, 1/2+ and 3/2+ , are not
removed but do increase slightly due to the imperfect COM factorization.
While the experimental excited states of the particle-attached nuclei 23 O and 23 F are not
well-known, the EOM states can be treated in an equivalent way to remove the negativeparity states in 23 O, shown in Fig. 4.10. The low-lying excited states of the particle-removed
nuclei 15 N and 15 O are shown in Fig. 4.11 and those of 21 N and 21 O are shown in Fig. 4.12.
These EOM-CC states can now be used to calculate properties of excited-states and openshell systems. In particular, the transition amplitudes between neutron- and proton-attached
states, or between neutron- and protron-removed states, describes the nuclear structure
components of the corresponding beta-decay.

99

17 F(5/2 + )

0.8

23 O(1/2 + )

2

23 O, 23 F

0.6

5

10 15 20 25 30
ω (MeV)

0.4
0.2

0

0.0
1.0
15 N(1/2 − )

15 N(1/2 − )

15 O(1/2 − )

15 O(1/2 − )

0.8

21 N(1/2 − )
21 O(5/2 + )

21 N(1/2 − )
21 O(5/2 + )

15 N, 15 O

21 N, 21 O

0.6

5

10 15 20 25 30
ω (MeV)

0.4

f

E cm (ω) (MeV)

6

18
17
16
15
14
13
12
11

f

10

E cm (ω) (MeV)

17 O, 17 F

f

4

8

23 O(1/2 + )
23 F(5/2 + )

E cm (ω) (MeV)

E cm (ω) (MeV)

6

23 F(5/2 + )

18
17
16
15
14
13
12

f

17 F(5/2 + )

ω (MeV)

8

17 O(5/2 + )

ω (MeV)

1.0
17 O(5/2 + )

4

0.2

2

0.0

0
10

15

20

ω (MeV)

25

30

10

15

20

ω (MeV)

25

30

Figure 4.8: Ground-state COM energies, Eq. (3.68), for open-shell nuclei at varies harmonic
oscillator frequencies with the NN+3N(400)-induced interaction with λSRG = 2.0 fm−1 at
emax = 12. The top row shows the results for the particle-attached nuclei 17 O, 17 F, 23 O,
and 23 F, and the bottom row shows the results for the particle-removed nuclei 15 N, 15 O,
21 N, and 21 O. The right column shows that the COM wave function practically vanishes
according to Eq. (3.67).

100

8

E (MeV)

6

5/2 +
1/2 +
3/2 +
3/2 −
1/2 −

5/2 +
1/2 +
3/2 +
3/2 −
1/2 −

17 O

17 F

4
2
0

β=0

β=1

β=0

Exp.

β=1

Exp.

Figure 4.9: Low-lying PA-EOM states for 17 O and 17 F with and without a Lawson-Gloeckner
term, along with the experimentally determined spectra. The negative-partity states are
COM contaminants and are removed by artificially raising the COM excitation energy with
the parameter β according to Eq. (3.72).

8
7
6

1/2 +
3/2 +
3/2 −
1/2 −

5/2 +
1/2 +
3/2 +
3/2 −

23 O

23 F

E (MeV)

5
4
3
2
1
0

β=0

β=1

β=0

Exp.

β=1

Exp.

Figure 4.10: Low-lying PA-EOM states for 23 O and 23 F with and without a LawsonGloeckner term. The negative-partiy states in 17 O are COM contaminants and are removed
by artificially raising the COM excitation energy with the Lawson-Gloeckner method, Eq.
(3.72). The excited states of these nuclei have not been experimentally detemined.

101

12
10

1/2 −
3/2 −
1/2 +

1/2 −
3/2 −
1/2 +

15 N

15 O

E (MeV)

8
6
4
2
0

β=0

β=1

β=0

Exp.

β=1

Exp.

Figure 4.11: Low-lying PR-EOM states for 15 N and 15 O with and without a lawson-gloeckner
term, and the experimentally determined spectra. The 1/2+ state is a COM contaminant
and is removed by artificially raising the COM excitation energy with the parameter β
according to Eq. (3.72).

10

E (MeV)

8

1/2 −
3/2 −
1/2 +

5/2 +
1/2 +
3/2 −
1/2 −

21 N

21 O

6
4
2
0

β=0

β=1

β=0

Exp.

β=1

Exp.

Figure 4.12: Low-lying PA-EOM states for 21 N and 21 O with and without a LawsonGloeckner term, Eq. (3.72), and the experimentally determined spectra. The negative-partiy
states in 17 O are COM contaminants and are removed by artificially raising the COM excitation energy. The excited states of these nuclei have not been experimentally detemined.

102

Chapter 5
Beta–Decay Effective Operators
Now with the ability to calculate open-shell states, properties of and transitions between
these states can be obtained. In particular, this chapter focuses on beta-decay transitions
between two particle-attached states and between two particle-removed states. However,
because these techniques are formulated in a general way, they can be applied to any type of
one-body operator and any type of EOM-CC state. First, this chapter describes beta-decay
processes in detail and the relevant operators are introduced. Then, these operators are
used to construct effective coupled cluster operators that take important correlations from
the ground-state wave function into account. Finally, the effective operators are applied to
EOM-CC states to calculate transition amplitudes for the corresponding processes. These
amplitudes can then be used to calculate observables like decay strengths and half-lives.

5.1

Beta-Decay Properties

Nuclear beta decay describes a class of radioactive decays of the atomic nucleus that result
as a consequence of the weak interaction. These processes involve the exchange of a W
boson which allows a quark within a proton or neutron to change type, thus converting
a neutron to a proton or vice versa. Additionally, this conversion is accompanied by an
electron-antineutrino or positron-neutrino pair which ensures the process conserves charge
and lepton number. This work focuses on the three most common types of weak processes
that occur within atomic nuclei: β − decay, β + decay, and electron capture (EC). The first,
103

β − decay, is the process whereby a down quark becomes an up quark, converting a neutron
to a proton along with the creation of an electron, or β − particle, and an antineutrino.
Schematically, this decay can be written as,

β − decay :

n −→ p + e− + ν̄e .

The corresponding mirror process is the β + decay, whereby an up quark becomes a down
quark, converting a proton to a neutron with the creation of a positron, or β + particle, and
a neutrino. This decay is represented by,

β + decay :

p −→ n + e+ + νe .

A closely related process is the electron capture, whereby an atomic electron interacts with
an up quark of a proton via a W boson, causing the conversion to a down quark and the
release of a neutrino. This process can be written as,

Electron capture :

p + e− −→ n + νe .

These three processes are schematically represented in Fig. 5.1.
The decay Q-value is a measure of the net energy of these processes. In free space,
only β − decay has a positive Q-value and can occur spontaneously. However, within a
nucleus, correlations between nucleons cause the relatively simple exchanges in Fig. 5.1 to
take on more complicated, higher-order processes involving pion exchanges between different
nucleons, see Fig. 5.2. This also has the effect of changing the energetics of the β + decay and
electron-capture processes such that their Q-values are positive and can occur within certain

104

n

p
u
u

d

n

−

νe

d

e

W−

u

d

+

νe

d

e

u

W+

d

W+

u

d

e−

u

p

n

νe

p

(a)

(b)

(c)

Figure 5.1: Schematic representations of the three free-space weak processes in this work:
β − decay (a), β + decay (b), and electron-capture (c). The coupling constant for the point
interaction vertex is the weak-interaction coupling constant gW .
ΨF

νe
p

A−1
π0 π−

e−

ΨF

νe
n

A−1
π0 π+

n
ΨI

n νe

A−1
π0

p
ΨI

(a)

e+

ΨF
π+
−
p e

ΨI

(b)

(c)

Figure 5.2: Schematic representations of a higher-order weak interactions involving pionexchange that occur within a nucleus, for β − decay (a), β + decay (b), and electron-capture
(c). These processes are not included in the impulse approximation. The coupling constant
for the point interaction vertex is the effective weak-interaction coupling constant GF .
nuclei. However, these complicated processes can be approximated by assuming that the
weak interaction occurs at very short length and time scales due to the large mass of the W
boson, mW . Known as the impulse approximation, this treatment of weak processes within
nuclei ignores any nucleon-nucleon interactions involved in weak decays, instead treating a
single nucleon as a spectator to the other nucleons during the weak interaction, shown in
Fig. 5.3. In addition, this motivates the simplification of the W boson exchanges in Fig.
5.1, with coupling constant gW , to a point interaction, shown in Figs. 5.3 and 5.2, with an
√ 2
effective coupling constant GF = 2gW
/8(mW c2 )2 [156].
These beta-decay processes can be further characterized by the character of their angular
momentum. Allowed beta decays involve an electron-antineutrino or positron-neutrino pair
105

ΨF

ΨF
p νe

n νe

e−

A−1

ΨF
e+

A−1
n

ΨI

νe

p

e−

A−1
p

ΨI

(a)

n

ΨI

(b)

(c)

Figure 5.3: Schematic representations of the impulse approximations to the different weak
processes within an A-body nucleus: β − decay (a), β + decay (b), and electron-capture (c).
The active nucleon doesn’t interact with the initial and final nuclei during the weak process.
The coupling constant for each interation vertex is the weak-interaction coupling constant
GF .
with no angular momentum, L = 0. In this case, there can be no parity change between
the initial and final nuclear states. Also, each pair of spin- 12 leptons carries a coupled spin
of either S = 0 or S = 1. Therefore, the initial and final nuclear states must have angular
momenta that only differ by 0 or 1, (∆J = 0, 1). Additionally, parity must be conserved
for these interactions. Fermi transitions (F) occur when the nuclear states are coupled to
leptons with S = 0, and Gamow-Teller transitions (GT) occur when the lepton have a
spin S = 1. The summary of these selection rules are shown in table 5.1. It should be
Decay Type

∆J = JF − JI

πF πI

Fermi
Gamow-Teller
Gamow-Teller

0
1(JF = 0 or JI = 0)
0, 1(JF > 0 and JI > 0)

+1
+1
+1

Table 5.1: Summary of the selection rules for allowed beta decays according to the angular
momentum J and parity π of the initial (I) and final (F ) states.

stated that while Gamow-Teller transitions can occur without a change in the nuclear spin
(∆J = 0 for JF > 0 and JI > 0), the Gamow-Teller operator carries an angular momentum
of J = 1.
The half-life of these processes T1/2 can be calculated using a combination of both the
106

Fermi and Gamow-Teller type transitions in the form of their reduced transition amplitudes,
BF and BGT , respectively,
gV2 |MF |2
BF =
2Ji + 1

2 |M
2
gA
GT |
BGT =
,
2Ji + 1

(5.1)

where Ji is the final state angular momentum. The factor gV is the vector coupling constant,
and its value can be shown to be exactly gV = 1.0. The axial-vector coupling constant gA
has a free-space value of gA = −gV , but is altered within nuclei due to nucleon-nucleon
correlations. The exact problem of how to treat the value of the axial-vector constant has
been a widely studied topic for decades [157, 158, 159, 160], but this work will use the value
gA /gV = 1.261(8) [161]. The transition matrix elements MF and MGT are measures of the
overlap integral between the initial and final nuclear states for the different transitions and
will be discussed in the next section. Inserting these reduced transition amplitudes into the
standard result from time-dependent perturbation theory gives the decay half-life,

T1/2 =

f
.
K0 (BF + BGT )

(5.2)

The factor f represents a phase-space integral over the final nuclear and lepton states and
depends on the decay Q-value. The factor K0 encodes the relevant constants involved,

K0 =

2π 3 ~7 ln 2
≈ 6147 s,
m5e c4 G2F

(5.3)

where me is the electron mass, and GF is the effective coupling constant. The next section
describes the process of improving upon the impulse approximation by using coupled cluster
theory to include higher-body interactions with the CC similarity transformation.

107

5.2

Coupled Cluster Effective Operators

The main step in calculating dynamic properties within an ab initio framework is to calculate
the transition matrix elements of the relevant operator between two correlated many-body
states. In the cases considered in this thesis, the correlated many-body states are given by
right and left expansions for the PA-EOM(2) and PR-EOM(2) operators in Eqs. (4.3) and
(4.3), respectively. Beta-decay properties are computed with the Fermi and Gamow-Teller
operators which are one-body operators that change a neutron to a proton or vice versa and
carry the quantum numbers that correspond to the rules in table 5.1.
In the impulse approximation, the Fermi operator, which has no spin component, is
simply equivalent to the isospin raising operator for the β − Fermi transition, which changes
a neutron to a proton, or the lowering operator for the β + transition, which changes a proton
to a neutron,
ÔF ∓ =

X
pq

n
o
†
hpk τ̂± kqi âp âq .

(5.4)

The reduced matrix element, see appendix D, hpk τ̂± kqi, is given by, see [162],
hpk τ̂± kqi =

p

2jp + 1 δnp nq δlp lq δjp jq δtp tq ¹1 ,

(5.5)

where the quantum numbers n and l can refer to the quantum numbers of any spherical
basis, such as the harmonic oscillator basis, and t is the states isosping projection. Similarly,
the Gamow-Teller operator also includes the isospin raising/lowering operator in addition to
the spin operator σ̂,
ÔGT ∓ =

X
pq

n
o
†
hpk σ̂τ̂± kqi âp âq .

108

(5.6)

Here, the reduced matrix element, hpk σ̂τ̂± kqi, is also given by, see [162],
hpk σ̂τ̂± kqi =

q

(
6(2jp + 1)(2jq + 1)

1
2
jp

1
2
jq

1
lp

)

3

(−1)lp +jp + 2 δnp nq δlp lq δtp tq ±1 .

(5.7)

These operators can now be used to calculate the transition matrix elements in Eq. (5.1)
by applying them between an initial state, |Ψi i, and a final state, hΨf |. The Fermi reduced
transition amplitude is,

MF = hΨf k ÔF ∓ kΨi i = δJ Ji
f

X
pq

n
o
†
hpk τ̂± kqi hΨf k âp âq kΨi i .

(5.8)

Similarly, the Gamow-Teller reduced transition amplitude is,

MGT = hΨf k ÔGT ∓ kΨi i =

X
pq

n
o
†
hpk σ̂τ̂± kqi hΨf k âp âq kΨi i .

(5.9)

Within the coupled cluster framework, the initial and final states take the form of PAEOM or PR-EOM states given generically in Eq. (4.1). According to Eqs. (4.7) and (4.12),
the left and right eigenstates of the bare Hamiltonian can be given by hΨf | = hΦ0 | L̂f e−T̂
and |Ψi i = eT̂ R̂i |Φ0 i. Inserting these EOM states into Eqs. (5.8) and (5.9), gives,
MF = δJ Ji
f

X
pq

n
o
†
hpk τ̂± kqi hΦ0 k L̂f e−T̂ âp âq eT̂ R̂i kΦ0 i .

(5.10)

Similarly, the Gamow-Teller matrix element becomes,

MGT =

X
pq

n
o
†
hpk σ̂τ̂± kqi hΦ0 k L̂f e−T̂ âp âq eT̂ R̂i kΦ0 i .

(5.11)

The resemblance of these equations to Eq. (3.6) motivates the construction of an effective

109

operator, λ Ō. The Fermi and Gamow-Teller effective operators have the form,

ŌF ∓ =

X

ŌGT ∓ =

X

pq

pq

hpk τ̂±

kqi e−T̂

hpk σ̂τ̂±

n
o
†
âp âq eT̂ ,

kqi e−T̂

n
o
†
âp âq eT̂ .

(5.12)
(5.13)

n
o
†
The similarity-transformed component, e−T̂ âp âq eT̂ , is known as the one-body density
matrix. Using the Baker-Campbell-Housedorf expansion like Eq. (3.7) and the reduction to
connected diagrams like Eq. (3.8), the one-body density matrix can be reduced to the form,

n
o
n
o 
†
†
e−T̂ âp âq eT̂ =
âp âq eT̂ .
c

(5.14)

The effective operators are best calculated using diagrammatic techniques like those used
for the effective Hamiltonian. In this case, the bare one-body operators are depicted by
while the effective one- and two-body operators are depicted by

the vertex type
and

, respectively.

The one-body beta-decay operators can be split into hp, hh, pp, and ph components
which are analogous to the one-body components of the effective Hamiltonian, see appendix
B. The hp component has no connected terms in Eq. (5.14), and so it is unchanged by the
similarity transformation,

=
i

a

i

ÎťOi
a

a

= Îť Oai .

110

(5.15)

The pp component is augmented by single excitations from the reference state. This component’s diagrammatic and algebraic expressions are,

a

a

a

=
b

+

k
b

b

ÎťOa
b

= λ Oba −

X

Îť O k ta .
b k

(5.16)

k

Similarly, the hh component is also augmented by single excitations from the reference state,

j

j

j

=
i

+

c

i

ÎťOi
j

= Îť Oji +

X

i

Îť O i tc .
c j

(5.17)

c

The ph component of the effective beta-decay operator includes effects from single and double
excitations from the reference state. The diagrammatic and algebraic expressions for this
component are,

a

i

a

a

i

=

+

i

i

+

a

c

ÎťOa
i

= Îť Oia +

X
c

Îť O a tc
c i

−

k

X

Îť O k ta
i k

k

a

+

+

X

c

Îť O k tac .
c ik

i

k

(5.18)

kc

Like the higher-body interactions induced by the CC similarity transformation, two-body
effective operators are induced from the one-body bare beta-decay operator. Calculating
properties using PA-EOM(2) and PR-EOM(2) states requires only two components of the

111

effective two-body operator. The pphp component is the results of double excitations from
the reference state and is given by the following diagram and the corresponding expression,

b

a

b

i

a

=
c

c

Îť O ab
ic

=−

X

i

k

Îť O k tab .
c ik

(5.19)

k

Similarly, the hphh component is represented by the following diagram and its corresponding
algebraic expression,

j

k

j

a

k

=
i

i

Îť O ia
jk

=

X

a

c

Îť O i tca .
c jk

(5.20)

c

After constructing the effective operators, the reduced transition amplitudes in Eq. (5.1)
can be calculated. Because of the ambiguity in the bi-orthonormalization discussed in section
4.2, the square-norm of the matrix elements, |M |2 , must be written using both the left and
right solutions for each of the initial and final states. Expanded in terms of the EOM states,
these reduced amplitudes for the Fermi and Gamow-Teller operators become,

BF =

X
gV2
δJ J i
hΦ0 k L̂i ŌF ± R̂f kΦ0 i hΦ0 k L̂f ŌF ∓ R̂i kΦ0 i ,
2Ji + 1 f

(5.21)

f

BGT =

2
gA

X

2Ji + 1

f

hΦ0 k L̂i ŌGT ± R̂f kΦ0 i hΦ0 k L̂f ŌGT ∓ R̂i kΦ0 i .

(5.22)

Using the machinery developed in this thesis and the techniques discussed in this chapter,

112

the calculated half-lives of various nuclei will be presented in an upcoming paper.

113

Chapter 6
Conclusions and Perspectives

6.1

Summary and Conclusions

Accurate ab initio calculations of beta-decay transition amplitudes are necessary for answering many open questions from a wide range of areas in modern physics from nucleosynthesis
to fundamental symmetries. The accuracy and scope needed to answer such questions necessitate a technique that is widely applicable, systematically improvable, and scalable to large
systems. In this thesis, we have developed the formalism for and achieved the application of
techniques based on coupled cluster theory that fulfill these requirements.
The large and versatile program that implements these techniques can be used for many
different fermionic systems including the homogeneous electron gas, quantum dots, and finite nuclei. Because of the program’s modular form, additional systems like neutron drops,
infinite nuclear matter, and atomic systems can easily be added in subsequent updates. Importantly, we extended the single-determinant coupled-cluster method to open-shell systems
using the equation-of-motion method which grants a broader reach across the nuclear chart.
Also, we added the crucial ability to perform calculations with and without three-body forces
which ensures accurate results. Also modular, these components of the program can be easily extended to higher-order EOM approximations, two-particle-attached and two-particleremoved EOM states, and the inclusion of full three-body forces. Lastly, we’ve implemented
the ability to construct any effective one-body operator and calculate the corresponding

114

observables using ground and excited EOM states. In future iterations of this code, higherorder effective operators, like those required for double-beta-decay experiments, can also be
implemented. Also, these higher-order operators can be constructed from two-body chiral
weak currents and used to investigate the quenching of the axial vector coupling constant.
In addition to developing and implementing these various techniques, we performed calculations at each step to verify the results. In particular, we provided a proof of principle by
comparing our results with those from other ab initio methods for various different systems.
Also, we calculated ground-state energies as well as particle-attached and particle-removed
spectra for various light nuclei, focusing mainly on the oxygen chain. In future publications,
this machinery will be extended up the nuclear chart to heavier nuclei and out to the limits
of stability to calculate beta-decay properties of nuclei around 78 Ni and 100 Sn, which are
important to future experiments at FRIB. For example, consistent beta-decay lifetime calculations from ab initio methods will be invaluable for astrophysical simulations of different
nucleosynthesis processes.

115

APPENDICES

116

Appendix A
CCSD Diagrams
The following diagrams and their corresponding algebraic expressions comprise the different contributions to the CCSD cluster amplitudes without directly building the effective
Hamiltonian, H . The boxed diagrams are automatically zero in a Hartree-Fock basis. The
contributions to the CCSD singles equation, (3.10), are given by Eqs. (A.1)-(A.8).

a

i

i

fˆN t̂1 |Φ0 ic =

= fia −

a

a

+

+
k

X

fik tak +

X

i
c

fca tci

(A.1)

c

k

a

V̂N t̂1 |Φ0 ic =

c

k

=−

X

117

kc

c
Vka
ic tk

i

(A.2)

fˆN t̂2 |Φ0 ic =

=

a
c

i

k

P k ac
fc tki

(A.3)

kc

a

V̂N t̂2 |Φ0 ic =

=

i
i

c

k

a

+

d

c

k

l

1 X ka cd 1 X kl ca
Vcd tki +
Vic tkl
2
2
kcd

(A.4)

klc

i

fˆN t̂21 |Φ0 ic

a

=

=

l

d

P l ad
fd tl ti

(A.5)

kcd

a

V̂N t̂21 |Φ0 ic =

=

i
i

c

X

k

d

c d
Vka
cd tk ti +

kcd

c

X
klc

118

a

+
k

c a
Vkl
ic tk tl

l

(A.6)

i

V̂N t̂1 t̂2 |Φ0 ic =

a

a

=−

l

i

+

d

c

d

k

a

+

l

c

k

c

k d

1 X kl cd a 1 X kl ca d X kl ad c
Vcd tki tl −
Vcd tkl ti +
Vcd til tk
2
2
klcd

klcd

i

l

(A.7)

klcd

i

V̂N t̂31 |Φ0 ic =

a
l

=−

d

X

c

k

c d a
Vkl
cd tk ti tl

(A.8)

klcd

The contributions to the CCSD doubles equation, (3.10), are given by Eqs. (A.9)-(A.19).

a

i

b

j

V̂N |Φ0 ic =

= Vab
ij

(A.9)

b

fˆN t̂2 |Φ0 ic =

i

a

= P̂ (ab)

j

j

+

a

i

b

c

X
c

fcb tac
ij − P̂ (ij)

119

k

X
k

fjk tab
ik

(A.10)

i

V̂N t̂1 |Φ0 ic =

b

V̂N t̂2 |Φ0 ic =

+

=

X

X

c
Vab
cj ti

b

i

+

(A.11)

c

l

i
j

c

+

b

a

j

d

k

c

X
1 X kl ab 1 X ab cd
ac
Vkb
Vij tkl +
Vcd tij − P̂ (ij|ab)
ic tkj
2
2

j

a

X
kl

a
b

+

cd

i
j

c

X

b

i

l

a b
Vkl
ij tk tl +

(A.12)

kc

cd

k

=

a
Vkb
ij tk + P̂ (ij)

a
b

i

=

b

k

kl

V̂N t̂21 |Φ0 ic

j

c

j

a
k

a
i

k

= − P̂ (ab)

i

j

a

+

d

c d
Vab
cd ti tj − P̂ (ij|ab)

120

b

a

j
k

X
kc

c

a c
Vkb
ic tk tj

(A.13)

a

V̂N t̂22 |Φ0 ic =

b

i
c

k

l

a

c

k

a

+

d

j

i
d

l

c

b

i

b

+

k

a

c

d

k

j

l

i
j

+

=

j

l

b

d

X
X
1 X kl ab cd
ac tbd − P̂ (ij) 1
Vcd tkl tij + P̂ (ab)
Vkl
t
Vkl
tab
tcd
cd
lj
ki
cd
lj
ki
4
2
klcd

klcd

klcd

1 X kl db ca
− P̂ (ab)
Vcd tij tkl
2

(A.14)

klcd

a

fˆN t̂1 t̂2 |Φ0 ic =

i

i
b

k

= − P̂ (ab)

c

j

+

a
j
c

P k a cb
P
fc tk tij − P̂ (ij) fck tci tab
kj
kc

kc

121

b

k

(A.15)

a

i

i
b

V̂N t̂1 t̂2 |Φ0 ic =

l

c

i

+

j

j

k

= P̂ (ij|ab)

+
d c

X
kcd

k

j

a
j

c

X
klc

b

+

k d

bc d
Vka
cd tjk ti − P̂ (ij|ab)

c

i
b

c

d
j

+

k

+

l

b

i

i

a
b

c

j

k

a

a

a

k

b

l

bc a
Vkl
ci tjk tl − P̂ (ab)

1 X kb cd a
Vcd tij tk
2
kcd

X
X
1 X kl ab c
tck tab
Vkl
− P̂ (ij)
+ P̂ (ij)
Vka
tck tdb
Vcj tkl ti + P̂ (ab)
ci
ij
lj
cd
2

a

V̂N t̂31 |Φ0 ic

=

i

i

d
j

k

c

= − P̂ (ij|ab)

+

b

X

a

j
b

c

k

l

a c d
Vkb
cd tk ti tj + P̂ (ij|ab)

kcd

X
klc

122

(A.16)

klc

kcd

klc

c a b
Vkl
cj ti tk tl

(A.17)

a

V̂N t̂21 t̂2 |Φ0 ic =

b

i
c

k

l

i

a

+

d

c

k

b

d

i
c

k

a
j

l

j

l

a

+

d

a
b

+

=

j

c

l

k

l

i
d

c

b
k

i
j

+

j

b

d

X
1 X kl ab c d 1 X kl cd a b
Vcd tkl ti tj +
Vcd tij tk tl + P̂ (ij|ab)
Vkl
tac
tbk tdi
cd
lj
2
2
klcd

klcd

klcd

1 X kl ab c d
1 X kl db c a
− P̂ (ij)
Vcd tlj tk ti − P̂ (ab)
Vcd tij tk tl
2
2
klcd

a

V̂N t̂41 |Φ0 ic

(A.18)

klcd

=

=

b
k

X

i
c

l

a b c d
Vkl
cd tk tl ti tj

klcd

123

j
d

(A.19)

Appendix B
Effective Hamiltonian Diagrams
The following diagrams and their corresponding algebraic expressions comprise the different
components of the CCSD effective Hamiltonian, H . Because some components are used as
intermediates to build other components, they must be built in the order written. Some
intermediate components overcount some diagrams which motivates the need for further
intermediates, denoted by X 0 , X 00 , and X 000 . The boxed diagrams are automatically zero in
a Hartree-Fock basis. The one-body components to the CCSD effective Hamiltonian, Eq.
(3.8), are given by Eqs. (B.1)-(B.4).

=
i

a

+
i

i

a

Xai = fai +

X

a

c

k

c
Vik
ac tk

(B.1)

kc

a

a

a

=
b

+
b

a

a
k

c

l

+

+
b

b

c

1 X kl ac X ka c X k a
Xba = fba −
Vbc tkl +
Vcb tk −
Xb tk
2
klc

kc

124

k

k

k
b

(B.2)

j

j

j

=
i

j

+
i

c

d

k

+
c

i

i

k

1 X ik cd X ik c
X 0ij = fji +
Vcd tjk +
Vjc tk
2
kcd

j

kc

j

j

=

+

i

(B.3)

c

i

i

Xji = X 0ij +

X

Xci tcj

(B.4)

c

Once the one-body components have been constructed, the pseudo-linear form for the
CCSD singles equation, Eq. (3.9), can be evaluated,

a

i

a

a

i

=0=

+

i

i

+

k

a

+

Xia = 0 = fia +

c

+

1X
2

kcd

c

Xca tci −

cd
Vka
cd tki −

a

+

d

X

c

k

i
i

k

i

+

c

c

a

a

X

k

c

a
X 0k
i tk +

X

ac
Vkl
ic tkl +

X

k

1X
2

l

a

+

i

k

c
Vka
ci tk

kc

klc

Xck tac
ik .

(B.5)

kc

During a CC iteration, it’s possible to use these updated singles amplitues T̂ 1 when evaluating
the doubles amplitudes T̂ 2 in the same iteration, which can accelerate the convergence.

125

The two-body components to the CCSD effective Hamiltonian, Eq. (3.8), are given by
Eqs. (B.6)-(B.20).

a

a

a

=
c

b

i

c

b

i

b

ia = Via −
Xbc
bc

X

a

i

0ij

ij

j

a

i

ij

ij

c

i

j

a

1 X ij c
Vca tk
2 c

(B.8)

k

+

a

(B.7)

1
+
2

k

j

i

k

=
i

b

k

X ka = Vka +

k

c

a
Vik
bc tk

k

j

k

i

=
i

i

(B.6)

+
c

k

b

a

a

b

c

k

=
c

k

1 X ik a
Vbc tk
2

ia
X 0ia
bc = Vbc −

a

i

1
+
2

j

Xka = Vka +

a

X ij
Vca tck
c

126

c

i

j

a

(B.9)

a

b

a

a

b

=
c

+

d

c

b

a

X

c

(B.10)

a

d

ab = X 0ab +
Xcd
cd

d

a
X 0kb
cd tk

+

d

c

k

b

=
c

k

d

ab
X 0ab
cd = Vcd − P̂ (ab)

a

b

b
k

l

c

d

1 X kl ab
Vcd tkl
2

(B.11)

kl

k

l

k

=
i

j

+
i

Xkl = Vkl +

l
c

j

ij

ij

k

l

l

k

+

d

i

i

j

j

c

X 0ij
1 X ij cd
Vcd tkl + P̂ (kl)
X kc tcl
2
c

(B.12)

cd

j

a

j

=
i

b

j

a

+
i

ia
X 0ia
jb = Vjb +

b

X
c

c
X 0ia
cb tj −

127

a
c

i

b

1 X ik a
Vjb tk
2
k

a

j

1
+
2i

b

k

(B.13)

j

a

j

a

i

b

j

=
i

b

ia
X 00ia
jb = Vjb +

a

1
+
2

c

i

a

j

1
+
2i

b

k

b

1 X 0ia c 1 X ik a
X cb tj −
Vjb tk
2 c
2

(B.14)

k

j

a

j

j

a

=
i

b

i

b

ia
X 000ia
jb = Vjb +

a

1
+
2

a

j

+

c

i

b

i

k

b

1 X ia c X ik a
Vjb tk
Xcb tj −
2 c

(B.15)

k

j

a

j

 j
1
+
2

a

=
i

b

i

ia
Xjb

The factor of

 
1
2

=

b

X 000ia
jb

a
c

j

1
+
2

k

i

a
c

i

b

b

 X
X
1
ia tc
ca + 1
Vik
Xcb
t
−
j
cb jk
2
2 c
kc

(B.16)

is applied when solving the CCSD equations but omitted when applying

the effective Hamiltonian in post-CC methods.

b

a

i

b

a

i

=
c

c

ab
X 0ab
ic = Vic +

b

1
+
2c

a

b
i
d

a

+
c

X
1 X ab d
a
Vdc ti − P̂ (ab)
X 00kb
ic tk
2
k

d

128

i

k

(B.17)

b

a

i

b

a

b

i

=
c

a

+

b

ab = Vab +
Xic
ic

b
a

+

X

d

X

k

a

j

X

j

k

a

k

=
i

j

a

a

k
i

k

+
c

k

a

+

j

X

c

a
Vil
jk tl + P̂ (jk)

d

i

lc

X

a

c

c
X 000ia
jc tk

c

l

X

k

+

l
i

ia = Via −
Xjk
jk

+ P̂ (jk)

a

a

i

j

c

(B.19)

j

l

i

+

a
l

l

k

j

k

1
+
2i

+
i

(B.18)

1 X il a
Vjk tl
2

j

a

1 X kl ab
Xic tkl
2
kl

i

j

k

a
X 0kb
ic tk

kb tad +
Xdc
ik

k

ia
X 0ia
jk = Vjk −

l

k

=
i

+

k

kd

k

i

a

i

c

d
Vab
dc ti − P̂ (ab)

Xck tab
ik + P̂ (ab)

j

k

c

c

d

X

i

k

c

a

b
a

+

i

+

d

c

c

−

b
i

il tca +
Xjc
lk

1X
2

cd

129

ia tcd +
Xcd
jk

X
c

Xci tca
jk

(B.20)

a

i

b

j

a

i

b

a

j

=0=

i
i

+

b

j

a

+

c

a

+

b

i

i
j

c

+

j

a

d

i
b

k

+

l

ab = 0 = Vab + P̂ (ab)
Xij
ij

X
c

+

b

k

a

b

a

j
k

j

+

c

Xca tcb
ij − P̂ (ij)

i

j

a

i

j

b

+
k

X

b

c

Xik tab
kj

k

X
X
X
1 X 0ab cd 1 X kl ab
kb tac − P̂ (ab)
c
X cd tij +
Xij tkl − P̂ (ab|ij)
Xic
X 0kb
tak + P̂ (ij)
X 0ab
ij
ic tj
kj
2
2
c
cd

kl

kc

k

(B.21)

130

Appendix C
Computational Implementation
The sums involved in building the CC effective Hamiltonian, solving the CC equations,
solving the EOM-CC equations, and building effective operators can all be reformulated
as matrix-matrix multiplications and thus performed with efficient LAPACK and BLAS
routines. To take advantage of this efficiency, the various cluster amplitudes and matrix
elements must be grouped into structures with similar index organization so that summed
indices map to the same states and matrix elements. An additional benefit to these structures
is that angular-momentum-coupling coefficients are automatically removed by summing over
Clebsch-Gordon coefficients, see chapter D.

Symmetry Channels
Each matrix structure is separated into different symmetry channels for different perutations
of its indices. The channels are denoted as ΣΞ~, where Ξ~ represents the relevant quantum
numbers of a certain channel. There are four different channel types that are relevant for
the structures used in this work.
The direct two-body channel categorizes the vector sum of two single-particle-state quantum numbers, ΣΞ~ ,
1
Ξ~pq = Ξ~p + Ξ~q

−→

|pqi ∈ ΣΞ~ =Ξ~ .
1 pq

(C.1)

The cross two-body channel categorizes the vector difference of two single-particle-state

131

quantum numbers or, equivalently, the vector sum of a the quantum numbers of a singleparticle state and a time-reversed single-particle state, ΣΞ~ ,
2
ξ~pq̄ = ξ~p − ξ~q = ξ~p + ξ~q̄

−→

|pq̄i ∈ Σξ~ .
2

(C.2)

The one-body channel categorizes single-particle states by their vector quantum numbers,
ΣΞ~ ,
3

ξ~p −→

|pi ∈ ΣΞ~ .

(C.3)

3

The cross three-body state categorizes the vector difference between the quantum numbers
of a direct two-body state and a single-particle state or, equivalently, the vector sum of the
quantum numbers of a two-body direct state and a time-reversed single-particle state, ΣΞ~ ,
3

ξ~pqr̄ = ξ~p + ξ~q − ξ~r = ξ~p + ξ~q + ξ~r̄ −→

|pqr̄i ∈ Σξ~ =ξ~ .
3 pqr̄

Channel-Partitioned Structures
Different matrix structures are indexed by their channel type: 1 for direct channels, 2 for cross
channels, and 3 for one-/three-body channels. For matrices with more than one structure
of the same type, there is an additional index that depends on the specific permutation
involved.
n
o
p
†
For a one-body operator Aq âp âq , there is a direct-channel matrix element and a

132

cross-channel matrix element,

p

A2 = Apq̄ .

A1 = Aq ,

(C.5)

n
o
pq
† †
For a two-body operator Ars âp âq âs âr , there is a direct-channel matrix element, four
cross-channel matrix elements, and four one-channel matrix elements,

pq

A1 = Ars ,
qr̄

ps̄

A22 = Asp̄ ,

A31 = Arsq̄ ,

A32 = Arsp̄ ,

A21 = Arq̄ ,
p

q

pq

For an EOM operator of the form Ar

pr̄

A23 = Asq̄ ,
pqs̄

A33 = Ar ,

qs̄

A24 = Arp̄ ,
pqr̄

A34 = As .

(C.6)

n
o
† †
âp âq âr , there is a direct-channel matrix ele-

ment, a one-channel matrix element, and two cross-channel matrix elements,

pq

A3 = Apqr̄ ,

p

A22 = Arp̄ .

A1 = Ar ,

q

A21 = Arq̄ ,

(C.7)

n
o
p
†
EOM operators of the form Aqr âp âr âq have similar structures,

p

A3 = Aqrp̄ ,

pr̄

A22 = Ar .

A1 = Aqr ,

pq̄

A21 = Aq ,

133

(C.8)



T̂
Matrix Form of H = ĤN e

c

tck̄
Xai = fai + Viā
ck̄
hp

hp

hhpp

X2 ←− f2 + V2
3

¡ t2

(C.9)

1
ab̄ ck̄
a k
Xba = fba − taklc̄ Vklc̄
b + Vck̄ t − tk X b
2
1
pp
hp
hhpp
X3 ←− − t3 · V3
− t3 · X 3
3
2 1
pp

pp

hppp

X2 ←− f2 + V2
4

¡ t2

(C.10)

1
k̄ + Vij̄ tck̄
X 0ij = fji + Vicdk̄ tcd
j
ck̄
2
1 hhpp
hh
X0 3 ←− V3
¡ t3
3
2 1
hh

hhhp

X0 2 ←− f2hh + V2
3

¡ t2

(C.11)

1
ij̄ ck̄
k̄
i c
Xji = fji + Vicdk̄ tcd
j + Vck̄ t + X c tj
2
1 hhpp
hp
Xhh
3 ←− 2 V31 · t33 + X 3 · t3
hh

hhhp
3

X0 2 ←− f2hh + V2

134

¡ t2

(C.12)

1 a cdk̄ 1 a klc̄
aī ck̄
aī kc̄
Xia = fia + Xca tci − tak X 0k
i − Vck̄ t + 2 Vcdk̄ ti − 2 tklc̄ Vi + tkc̄ X
1 hppp
1
hh
ph
pp
hhhp
X3 ←− X3 · t3 − t3 · X0 3 + V3
· t3 − t3 · V3
4
3
2 2
2 1
ph

ph

hphp

X2 ←− f2 − V2
2

hp

¡ t2 + t2 ¡ X 2

(C.13)

3

1 a k
ia
X 0ia
bc = Vbc − 2 tk Vbcī
1
hppp
hppp
hhpp
− t3 · V 3
X0 3 ←− V3
2
2
2
2

(C.14)

ia = Via − ta Vk
Xbc
bc
k bcī
hppp
2

X3

hppp
2

←− V3

hhpp
2

− t3 · V3

1 ijā
0ij
ij
X ka = Vka + Vc tck
2
1 hhpp
hhhp
ij
¡ t3
X0 3
←− Vka + V3
3
2 3

ij

ij

(C.15)

(C.16)

ijā

Xka = Vka + Vc tck
hhhp
3

X3

ij

hhpp
3

←− Vka + V3

135

¡ t3

(C.17)

ab
a 0k
X 0ab
cd = Vcd − P̂ (ab) tk X cdb̄
pppp

X0 1

pppp
1(2)

X0 3

pppp

←− V1

hppp
1

←− ∓t3 · X0 3

ab = X 0ab + 1 tab Vkl
Xcd
cd
2 kl cd
1
pppp
hhpp
pppp
X1
←− X0 1 + t1 · V1
2

(C.18)

(C.19)

1 ij
ij
ij
0ij k̄
Xkl = Vkl + Vcd tcd
+ P̂ (kl) X c tcl
kl
2
1 hhpp
¡ t1
Xhhhh
←− Vhhhh
+ V1
1
1
2
hhhp
4

Xhhhh
←− ∓X0 3
3
3(4)

¡ t3

(C.20)

1 a k
ia
0iab̄ c
X 0ia
jb = Vjb + X c tj − 2 tk Vjbī
hphp
1

←− V2

hphp
3

←− X0 3

hphp
2

1
hhhp
←− − t3 · V3
2
2

X0 2
X0 3
X0 3

hphp
1
hppp
3

136

¡ t3
(C.21)

1 0iab̄ c 1 a k
ia
X 00ia
jb = Vjb + 2 X c tj − 2 tk Vjbī
hphp

hphp

X00 2 ←− V2
1
1

1 hppp
hphp
X00 3 ←− X0 3 · t3
3
3
2
1
hphp
hhhp
X00 3 ←− − t3 · V3
2
2
2

(C.22)

1 iab̄ c
a k
ia
X 000ia
jb = Vjb + 2 X c tj − tk Vjbī
hphp

hphp

X000 2 ←− V2
1
1

1 hppp
hphp
X000 3 ←− X0 3 · t3
3
3
2
hphp

hhhp

X000 3 ←− −t3 · V3
2
2

ia
Xjb
hphp
X2
1

=

Via
jb

←−

+ Xciab̄ tcj

hphp
V2
1
hppp
3

←− X3

hphp
2

←− −t3 · V3

X3

 
1
−
Vicb̄k̄ tcjāk̄
2

 
1
hhpp
−
V 2 t2
1
1
2

hphp
3

X3

− tak Vkjbī

(C.23)

¡ t3
hhhp
2

(C.24)

1 abc̄ d
ab
a 00k
X 0ab
ic = Vic + 2 Vd ti − P̂ (ab) tk X icb̄
1 pppp
pphp
pphp
X0 3 ←− V3
+ V3
¡ t3
3
3
2 3
pphp
1(2)

X0 3

hphp

←− ∓t3 · X00 3
1

137

(C.25)

ab = Vab + Vabc̄ td − P̂ (ab) ta X 0k − tabī X k + P̂ (ab) taī X k d¯ + 1 tab X kl
Xic
c
ic
k icb̄
k
d i
k d¯ cb̄
2 kl ic
pphp
3

←− V3

pphp
1(2)

←− ∓t3 · X0 3

pphp
4

←− −t3 · X3

pphp
2(3)

←− ∓t2 · X2

pphp

1
hhhp
←− t1 · X1
2

X3

pppp
3

+ V3

¡ t3

hphp
1

X3

X3
X2

X1

pphp
3

hp

4

hppp
3

3

(C.26)

1 a l
ia
X 0ia
jk = Vjk − 2 tl Vjk ī
1
hphh
hphh
− t3 · Vhhhh
X0 3
←− V3
32
2
2
2

000iaj̄ c
tk

ia = Via − ta Vl + P̂ (jk) X
Xjk
c
jk
l jk ī
hphh
2

hphh
2

←− V3

hphh
3(4)

←− ∓X000 3 · t3
4

hphh
1(3)

←− ∓X2

hphh

1 hppp
←− X1
¡ t1
2

hphh
1

←− X3 · t3

X3
X3

X2

X1

X3

(C.27)

1 ia cd
ij̄ ¯l
tjk + Xci tcjkā
+ P̂ (jk) X ¯tckā
+ Xcd
cl
2

− t3 · Vhhhh
3
2

hphp

hhhp
3

¡ t2

3

hp

(C.28)

1

138

ab = Vab + P̂ (ab) X a tc − P̂ (ij) tabj̄ X k + 1 X 0ab tcd + 1 tab X kl
Xij
c ij b̄
ij
i
k
2 cd ij
2 kl ij
aj̄

ī c
− P̂ (ab|ij) tkc̄ Xikc̄
− P̂ (ab) tak X 0k
+ P̂ (ij) X 0ab
c tj
b̄
ij b̄

←− V1

pphh
1(2)

←− ±X3 · t3 ∓ t3 X0 3
1
1

pphh
3(4)

0
←− ∓t3 · Xhh
3 ÂąX 3

pphh
1(2)

←− −t2 · X2

pphh
3(4)

←− t2 · X2

X3
X3

X2
X2

pphh

1 pppp
1
+ X0 1 ¡ t1 + t1 ¡ Xhhhh
1
2
2

pphh

X1

hphh

pp

pphp
t
4 3

3

hphp
1

1

hphp
1

1

Matrix Form of



(C.29)

H NR̂A+1
Âľ


c

= ωµR̂A+1
Âľ

a X kc̄ − 1 X a r cdk̄
ωk ra = Xca rc + rkc̄
2 cdk̄
1 hppp
hp
pp
¡ r3
ωk r ←− X3 · r + r2 · X2 − X3
1
2 2

(C.30)

c − r ab X k + 1 X ab r cd − P̂ (ab) r b X kc̄ − 1 tabī Vk r cd¯l
ωk riab = −Xcabī rc + P̂ (ab) Xcb riā
i
k
kc̄ iā
2 cd i
2 k cdÂŻl
1
hhpp
pphp
ωk r3 ←− −X3
· r − t3 · V3
¡ r3
1
4
2 3

ωk r2

1(2)

pp

hphp

←− ∓X3 · r2 ± r2 · X2
2
2
1

1 pppp
ωk r1 ←− −r1 · Xhh
¡ r1
3 + 2 X1

(C.31)

139

Matrix Form of L̂A+1
H N = EµL̂A+1
Âľ
Âľ
One main difference for the left eigenproblem is that the disconnected term is computed as
an outer product rather than with matrix-matrix multiplication.

1
Ek la = lc Xac − lcdk̄ Xacdk̄
2
1
pp
pphp
Ek l ←− l · X3 − l3 · X3
4
2

(C.32)

i = P̂ (ab) l X ib̄ − l X c − X i lk + P̂ (ab) liā X c + 1 li X cd
Ek lab
a
c abī
c
k ab
b
2 cd ab
iā lkc̄ − 1 l tcd¯l Vk
− P̂ (ab) Xkc̄
b
2 cd¯l k abī
1
hhpp
hppp
Ek l3 ←− −l · X3
− l3 · t3 · V3
3
1
2
2

Ek l2

1(2)

pp

hp

hphp

←− ±l ⊗ X2 ∓ l2 · X3 ± X2
2
1

¡ l2

2

1
pppp
Ek l1 ←− −Xhh
3 ¡ l1 + 2 l1 ¡ X 1

Matrix Form of



H NR̂µA−1


c

(C.33)

= ωµR̂µA−1

1
ωk ri = −rk Xik + Xck̄ rick̄ − rklc̄ Xiklc̄
2
1
hp
hhhp
ωk r ←− −r · Xhh
3 + X 20 · r21 − 2 r3 · X 33

140

(C.34)

a = −r X k − P̂ (ij) r aī X k + X a r c + 1 r a X kl − P̂ (ij) X aī r ck̄ − 1 r Vkld¯tc
ωk rij
c ij
k j
k ijā
ck̄ j
2 kl ij
2 kld¯ c ijā
1
hphh
hhpp
ωk r3 ←− −r · X3
− r3 · V3
¡ t3
1
1
3
2

ωk r2

1(2)

hphp
2

←− ±r2 · Xhh
3 Âą X2
2

¡ r2

2

1
pp
ωk r1 ←− X3 · r1 + r1 · Xhhhh
1
2

(C.35)

Matrix Form of L̂A−1
H N = EµL̂µA−1
Âľ
Again, the disconnected term is computed as an outer product rather than with matrixmatrix multiplication.

1 i klc̄
Ek li = −Xki lk − Xklc̄
l
2
1 hphh
¡ l3
Ek l ←− −XHH
· l − X3
3
2 1

(C.36)

1 ij
ij
ij
ijā
j
Ek la = P̂ (ij) li Xaj̄ − Xk lk + lc Xac − P̂ (ij) Xk lakī + Xkl lakl
2
1 ijā
ÂŻ
j
− P̂ (ij) l Xack̄ī − Vc tckld¯lkld
ck̄
2
1 hhpp
hhhp
Ek l3 ←− −X3
· l − V3
¡ t3 ¡ l3
1
3
2 3
Ek l2

1(2)

hp

hphp

←− ±l ⊗ X 0 ± Xhh
3 ¡ l22 ¹ l22 ¡ X 22
2

1
pp
Ek l1 ←− l1 · X3 + Xhhhh
¡ l1
2 1

141

(C.37)

Appendix D
Angular Momentum Coupling
Before deriving useful equations for J-scheme angular momentum coupling, it’s necessary to
list some shorthand notations, definitions, and useful relationships:

p̂ ≡
X
{m}

p

2jp + 1

≡ sum over all m

(D.1)
(D.2)

Clebsch-Gordan coefficients:

hpq|Ji ≡ hjp mp ; jq mq |JM i

(D.3)

hpq̄|Ji ≡ hjp , mp ; jq , −mq |JM i(−1)(q−mq )

(D.4)

X
JM

hpq|Jihp0 q 0 |Ji = δmp m 0 δmq m 0
p
q

X
mp mq

hpq|Jihpq|J 0 i = δJJ 0 δM M 0

(D.5)
(D.6)

(D.7)

142

Clebsch-Gordan coefficient symmetries:

hjp mp ; jq mq |JM i = (−1)jp +jq −J hjp mp ; jq mq |JM i

(D.8)

= (−1)jp +jq −J hjq mq ; jp mp |JM i

(D.9)

Jˆ
= (−1)jp −mp hjp mp ; J − M |jq − mq i
q̂
Jˆ
= (−1)jq +mq hJ − M ; jq mq |jp − mp i
p̂
Jˆ
= (−1)jp −mp hJM ; jp − mp |jq mq i
q̂
Jˆ
= (−1)jq +mq hjq − mq ; JM |jp mp i
p̂

(D.10)
(D.11)
(D.12)
(D.13)

Six-J symbols:
(
p
s
(
X
j
ĵ32 1
j4

)
(
jp jq jr
≡
js jt ju
)(
)
δj j 0
j2 j3
j1 j2 j3
6
= 26
0
j5 j6
j4 j5 j6
ĵ3
j3
(
)
X
X p q J
hps̄|J 0 ihrq̄|J 0 i = Jˆ02
hpq|Jihrs|Ji
r s J0
0
JM
M
(
)
X
X p q J
0
0
hpq|Jihrs|Ji = Jˆ2
0 hps̄|J ihrq̄|J i
r
s
J
0 0
M

q r
t u

)

(D.14)

(D.15)

(D.16)

(D.17)

J M

Two-body, scalar J-scheme matrix elements (T̂ , Ĥ, H ), in terms of M-scheme matrix

143

elements:

pq J

X J =
rs
0
ps̄J
X 0
rq̄ J

=

p

pq J s̄

{m}

X
{m}

X J =
rs q̄
Xr

X

=

X
{m}

X
{m}

pmp qmq

Xrmr sms hpq|Jihrs|Ji
pmp qmq

Xrmr sms hps̄|J 0 ihrq̄|J 0 i
pmp qmq

Xrmr sms hrs|JihJ q̄|pi
pmp qmq

Xrmr sms hpq|JihJ s̄|ri

(D.18)

(D.19)

(D.20)

(D.21)

Two-body M-scheme matrix elements in terms of J-scheme, scalar matrix elements:

pmp qmq

Xrmr sms =

X
JM

=

pq J

X J hpq|Jihrs|Ji
rs

X
J 0M 0

=

X

=

X

JM

0
ps̄J
X 0 hps̄|J 0 ihrq̄|J 0 i
rq̄ J
p

X J hrs|JihJ q̄|pi
rs q̄
pq J s̄

Xr

JM

hpq|JihJ s̄|ri

(D.22)
(D.23)

(D.24)
(D.25)

To find the relationship between the scalar matrix elements of T̂ , Ĥ, and H with different
couplings, the M-scheme expressions are written in terms of their different couplings, then the
Clebsch-Gordon coefficients are reorganized using Eqs. (D.8)–(D.17) so that they have the
same form. A few examples of this recoupling are shown below with the relevant equation
pmp qmq

used at each step. The shorthand X ≡ Xrmr sms is used for clarity. The relationship

144

0
ps̄J

pq J

between X J and X 0 is,
rs
rq̄ J
pq J
X J
rs

)
(
p
q
J
0
0
=
Xhpq|Jihrs|Ji =
X Jˆ2
0 hps̄|J ihrq̄|J i
r
s
J
M {m}
J 0 M 0 {m}
)
(
X ps̄J 0
p q J
(D.19)
=
X 0 Jˆ2
J
r s J0
rq̄
X

X

(D.17)

(D.26)

J

pq J

p

As another example, the relationship between X J and X J is,
rs
rs q̄
pq J

X J =
rs
=

X
{m}

X
{m}

Xhpq|Jihrs|Ji =

X

Xhqp|Jihrs|Ji(−1)jp +jq −J

{m}

Jˆ
X hJ q̄|pihrs|Ji(−1)jp +jq −J
p̂

Jˆ
p
= X J (−1)jp +jq −J
rs q̄ p̂

(D.9)

(D.19)

(D.12)

(D.27)

When sums are formulated in the terms of channel-partitioned matrices like those in
sections 3.4.2 and C, the factors related to a structure’s angluar momentum coupling are
automatically summed with the identity Eq. (D.6). To demonstrate this, an example is
shown here. First, from the CCSD equations the sum in Eq. (3.44) is rewritten in terms
of the J-scheme structures. The indices represent single-particle states in the M-scheme

145

expression but represent degenerate shells in J-scheme,

J2 c̄
1 X kl db ca
1X a
Vcd tij tkl −→
t J hkl|J1 ihJ1 c̄|aiVkl
hkl|J2 ihJ2 c̄|ditd J hij|J3 ihJ3 b̄|di
d
1
kl c̄
ij 3 b̄
2
2
klcd

=

klcd
J1 J2 J3
{m}





1X

X
J2 c̄ d

ta J Vkl
t
hkl|J1 ihkl|J2 i hJ1 c̄|aihJ2 c̄|dihij|J3 ihJ3 b̄|di
kl 1 c̄ d
ij J3 b̄
2
mk ml
klcd
J1 J2 J3
{m}
=

h
i
J2 c̄ d
1X a
t
δ
t J Vkl
δ
J1 J2 M1 M2 hJ1 c̄|aihJ2 c̄|dihij|J3 ihJ3 b̄|di
kl 1 c̄ d
ij J3 b̄
2
klcd
J1 J2 J3
{m}

=



1X
2

klcd
J1 J3
{m}

=

J1 c̄ d
ta J Vkl
t J 
d
kl 1 c̄
ij 3 b̄


X

M1 mc

hJ1 c̄|aihJ1 c̄|di hij|J3 ihJ3 b̄|di

h
i
J1 c̄ d
1X a
t J Vkl
t
δ
δ
ja jd ma md hij|J3 ihJ3 b̄|di
kl 1 c̄ d
ij J3 b̄
2
klcd
J1 J3
{m}

=

J1 c̄ d
1X a
t J Vkl
t J hij|J3 ihJ3 b̄|ai
kl 1 c̄ d
ij 3 b̄
2

(D.28)

klcd
J1 J3
{m}

This final form has the same structure as Eq. (D.24) so that the sum can be collected into
the following structure with no angular-momenutm coupling coefficients,

ta J

ij b̄

←

J1 c̄ d
1X a
t J .
t J Vkl
ij 3 b̄
kl 1 c̄ d
2

(D.29)

klcd
J1 J3

Amplitudes of different coupling structures are then gathered using relationships like those
in Eq. (D.26) and (D.27).

146

REFERENCES

147

REFERENCES
[1] H. Hergert, “Private communication.”
[2] M. Shao, H. M. Aktulga, C. Yang, E. G. Ng, P. Maris, and J. P. Vary, “Accelerating
nuclear configuration interaction calculations through a preconditioned block iterative
eigensolver,” CoRR, vol. abs/1609.01689, 2016.
[3] P. López Rı́os, A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, “Inhomogeneous backflow transformations in quantum monte carlo calculations,” Phys. Rev. E,
vol. 74, p. 066701, Dec 2006.
[4] R. Machleidt, “Chiral symmetry and the nucleon-nucleon interaction,” Symmetry,
vol. 8, no. 4, 2016.
[5] K. A. Brueckner and C. A. Levinson, “Approximate reduction of the many-body problem for strongly interacting particles to a problem of self-consistent fields,” Phys. Rev.,
vol. 97, pp. 1344–1352, Mar 1955.
[6] H. A. Bethe, “Nuclear many-body problem,” Phys. Rev., vol. 103, pp. 1353–1390, Sep
1956.
[7] J. Goldstone, “Derivation of the brueckner many-body theory,” Proceedings of the
Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 239,
no. 1217, pp. 267–279, 1957.
[8] F. Coester, “Bound states of a many-particle system,” Nuclear Physics, vol. 7, no. Supplement C, pp. 421 – 424, 1958.
[9] F. Coester and H. Kmmel, “Short-range correlations in nuclear wave functions,” Nuclear Physics, vol. 17, no. Supplement C, pp. 477 – 485, 1960.
[10] H. Kümmel, K. Lührmann, and J. Zabolitzky, “Many-fermion theory in exps- (or
coupled cluster) form,” Physics Reports, vol. 36, no. 1, pp. 1 – 63, 1978.
[11] J. C. Slater, “The theory of complex spectra,” Phys. Rev., vol. 34, pp. 1293–1322, Nov
1929.
[12] E. U. Condon, “The theory of complex spectra,” Phys. Rev., vol. 36, pp. 1121–1133,
Oct 1930.
[13] R. F. Bacher, “The interaction of configurations: sd−p2 ,” Phys. Rev., vol. 43, pp. 264–
269, Feb 1933.
[14] C. W. Ufford, “Configuration interaction in complex spectra,” Phys. Rev., vol. 44,
pp. 732–739, Nov 1933.

148

[15] J. Hubbard, “The description of collective motions in terms of many-body perturbation theory,” Proceedings of the Royal Society of London. Series A, Mathematical and
Physical Sciences, vol. 240, no. 1223, pp. 539–560, 1957.
[16] N. Hugenholtz, “Perturbation theory of large quantum systems,” Physica, vol. 23,
no. 1, pp. 481 – 532, 1957.
[17] H. Schaefer, Quantum chemistry: the development of ab initio methods in molecular
electronic structure theory. Oxford science publications, Clarendon Press, 1984.
[18] I. Shavitt and R. J. Bartlett, Many-Body Methods in Chemistry and Physics: MBPT
and Coupled-Cluster Theory. Cambridge Molecular Science, Cambridge University
Press, 2009.
[19] J. Čı́žek, “On the correlation problem in atomic and molecular systems. calculation
of wavefunction components in urselltype expansion using quantumfield theoretical
methods,” The Journal of Chemical Physics, vol. 45, no. 11, pp. 4256–4266, 1966.
[20] J. Čı́žek and J. Paldus, “Correlation problems in atomic and molecular systems iii. rederivation of the coupled-pair many-electron theory using the traditional quantum chemical methodst,” International Journal of Quantum Chemistry, vol. 5, no. 4, pp. 359–
379, 1971.
[21] J. Čı́žek and J. Paldus, “Coupled cluster approach,” Physica Scripta, vol. 21, no. 3-4,
p. 251, 1980.
[22] P. Piecuch, K. Kowalski, I. S. O. Pimienta, and M. J. Mcguire, “Recent advances
in electronic structure theory: Method of moments of coupled-cluster equations and
renormalized coupled-cluster approaches,” International Reviews in Physical Chemistry, vol. 21, no. 4, pp. 527–655, 2002.
[23] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, “Modern theory of nuclear forces,”
Rev. Mod. Phys., vol. 81, pp. 1773–1825, Dec 2009.
[24] R. Machleidt and D. Entem, “Chiral effective field theory and nuclear forces,” Physics
Reports, vol. 503, no. 1, pp. 1 – 75, 2011.
[25] S. Bogner, R. Furnstahl, and A. Schwenk, “From low-momentum interactions to nuclear structure,” Progress in Particle and Nuclear Physics, vol. 65, no. 1, pp. 94 – 147,
2010.
[26] R. Roth, J. Langhammer, A. Calci, S. Binder, and P. Navrátil, “Similarity-transformed
chiral N N + 3N interactions for the ab initio description of 12 C and 16 O,” Phys. Rev.
Lett., vol. 107, p. 072501, Aug 2011.
[27] P. Navrátil, J. P. Vary, and B. R. Barrett, “Large-basis ab initio no-core shell model
and its application to 12 C,” Phys. Rev. C, vol. 62, p. 054311, Oct 2000.

149

[28] P. Navrátil, S. Quaglioni, I. Stetcu, and B. R. Barrett, “Recent developments in nocore shell-model calculations,” Journal of Physics G: Nuclear and Particle Physics,
vol. 36, no. 8, p. 083101, 2009.
[29] B. R. Barrett, P. Navrátil, and J. P. Vary, “Ab initio no core shell model,” Progress in
Particle and Nuclear Physics, vol. 69, no. Supplement C, pp. 131 – 181, 2013.
[30] B. S. Pudliner, V. R. Pandharipande, J. Carlson, S. C. Pieper, and R. B. Wiringa,
“Quantum monte carlo calculations of nuclei with a ≤ 7,” Phys. Rev. C, vol. 56,
pp. 1720–1750, Oct 1997.
[31] S. C. Pieper, , and R. B. Wiringa, “Quantum monte carlo calculations of light nuclei,”
Annual Review of Nuclear and Particle Science, vol. 51, no. 1, pp. 53–90, 2001.
[32] J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt, and
R. B. Wiringa, “Quantum monte carlo methods for nuclear physics,” Rev. Mod. Phys.,
vol. 87, pp. 1067–1118, Sep 2015.
[33] K. Tsukiyama, S. K. Bogner, and A. Schwenk, “In-medium similarity renormalization
group for nuclei,” Phys. Rev. Lett., vol. 106, p. 222502, Jun 2011.
[34] K. Tsukiyama, S. K. Bogner, and A. Schwenk, “In-medium similarity renormalization
group for open-shell nuclei,” Phys. Rev. C, vol. 85, p. 061304, Jun 2012.
[35] H. Hergert, S. Binder, A. Calci, J. Langhammer, and R. Roth, “Ab initio calculations
of even oxygen isotopes with chiral two-plus-three-nucleon interactions,” Phys. Rev.
Lett., vol. 110, p. 242501, Jun 2013.
[36] S. K. Bogner, H. Hergert, J. D. Holt, A. Schwenk, S. Binder, A. Calci, J. Langhammer,
and R. Roth, “Nonperturbative shell-model interactions from the in-medium similarity
renormalization group,” Phys. Rev. Lett., vol. 113, p. 142501, Oct 2014.
[37] H. Hergert, S. K. Bogner, T. D. Morris, S. Binder, A. Calci, J. Langhammer, and
R. Roth, “Ab initio multireference in-medium similarity renormalization group calculations of even calcium and nickel isotopes,” Phys. Rev. C, vol. 90, p. 041302, Oct
2014.
[38] H. Hergert, “In-medium similarity renormalization group for closed and open-shell
nuclei,” Physica Scripta, vol. 92, no. 2, p. 023002, 2017.
[39] S. R. Stroberg, H. Hergert, J. D. Holt, S. K. Bogner, and A. Schwenk, “Ground and
excited states of doubly open-shell nuclei from ab initio valence-space hamiltonians,”
Phys. Rev. C, vol. 93, p. 051301, May 2016.
[40] S. R. Stroberg, A. Calci, H. Hergert, J. D. Holt, S. K. Bogner, R. Roth, and A. Schwenk,
“Nucleus-dependent valence-space approach to nuclear structure,” Phys. Rev. Lett.,
vol. 118, p. 032502, Jan 2017.

150

[41] V. Somà, C. Barbieri, and T. Duguet, “Ab initio gorkov-green’s function calculations
of open-shell nuclei,” Phys. Rev. C, vol. 87, p. 011303, Jan 2013.
[42] V. Somà, C. Barbieri, and T. Duguet, “Ab initio self-consistent gorkov-green’s function
calculations of semi-magic nuclei: Numerical implementation at second order with a
two-nucleon interaction,” Phys. Rev. C, vol. 89, p. 024323, Feb 2014.
[43] V. Somà, A. Cipollone, C. Barbieri, P. Navrátil, and T. Duguet, “Chiral two- and threenucleon forces along medium-mass isotope chains,” Phys. Rev. C, vol. 89, p. 061301,
Jun 2014.
[44] M. Wloch, D. J. Dean, J. R. Gour, M. Hjorth-Jensen, K. Kowalski, T. Papenbrock,
and P. Piecuch, “Ab-initio coupled-cluster study of 16 O,” Phys. Rev. Lett., vol. 94,
p. 212501, Jun 2005.
[45] M. Woch, J. R. Gour, P. Piecuch, D. J. Dean, M. Hjorth-Jensen, and T. Papenbrock,
“Coupled-cluster calculations for ground and excited states of closed- and open-shell
nuclei using methods of quantum chemistry,” Journal of Physics G: Nuclear and Particle Physics, vol. 31, no. 8, p. S1291, 2005.
[46] G. R. Jansen, J. Engel, G. Hagen, P. Navratil, and A. Signoracci, “Ab-initio coupledcluster effective interactions for the shell model: Application to neutron-rich oxygen
and carbon isotopes,” Phys. Rev. Lett., vol. 113, p. 142502, Oct 2014.
[47] G. R. Jansen, M. D. Schuster, A. Signoracci, G. Hagen, and P. Navrátil, “Open sd-shell
nuclei from first principles,” Phys. Rev. C, vol. 94, p. 011301, Jul 2016.
[48] G. Hagen, A. Ekström, C. Forssén, G. R. Jansen, W. Nazarewicz, T. Papenbrock, K. A.
Wendt, S. Bacca, N. Barnea, B. Carlsson, C. Drischler, K. Hebeler, M. Hjorth-Jensen,
M. Miorelli, G. Orlandini, A. Schwenk, and J. Simonis, “Neutron and weak-charge
distributions of the 48 Ca nucleus,” Nature Physics, vol. 12, pp. 186 EP –, Nov 2015.
Article.
[49] K. Kowalski, D. J. Dean, M. Hjorth-Jensen, T. Papenbrock, and P. Piecuch, “Coupled
cluster calculations of ground and excited states of nuclei,” Phys. Rev. Lett., vol. 92,
p. 132501, Apr 2004.
[50] J. R. Gour, P. Piecuch, M. Hjorth-Jensen, M. Wloch, and D. J. Dean, “Coupled-cluster
calculations for valence systems around 16 O,” Phys. Rev. C, vol. 74, p. 024310, Aug
2006.
[51] S. Binder, P. Piecuch, A. Calci, J. Langhammer, P. Navrátil, and R. Roth, “Extension
of coupled-cluster theory with a noniterative treatment of connected triply excited
clusters to three-body hamiltonians,” Phys. Rev. C, vol. 88, p. 054319, Nov 2013.
[52] T. D. Morris, J. Simonis, S. R. Stroberg, C. Stumpf, G. Hagen, J. D. Holt, G. R. Jansen,
T. Papenbrock, R. Roth, and A. Schwenk, “Structure of the lightest tin isotopes,” 2017.

151

[53] N. Cabibbo, “Unitary symmetry and leptonic decays,” Phys. Rev. Lett., vol. 10,
pp. 531–533, Jun 1963.
[54] M. Kobayashi and T. Maskawa, “Cp-violation in the renormalizable theory of weak
interaction,” Progress of Theoretical Physics, vol. 49, no. 2, pp. 652–657, 1973.
[55] I. S. Towner and J. C. Hardy, “The evaluation of v ud , experiment and theory,”
Journal of Physics G: Nuclear and Particle Physics, vol. 29, no. 1, p. 197, 2003.
[56] I. S. Towner and J. C. Hardy, “An Improved calculation of the isospin-symmetrybreaking corrections to superallowed Fermi beta decay,” Phys. Rev., vol. C77,
p. 025501, 2008.
[57] I. Towner, “Quenching of spin operators in the calculation of radiative corrections for
nuclear beta decay,” Physics Letters B, vol. 333, no. 1, pp. 13 – 16, 1994.
[58] I. Towner, “The nuclear-structure dependence of radiative corrections in superallowed
fermi beta-decay,” Nuclear Physics A, vol. 540, no. 3, pp. 478 – 500, 1992.
[59] F. Barker, B. Brown, W. Jaus, and G. Rasche, “Determination of —vud— from fermi
decays and the unitarity of the km-mixing matrix,” Nuclear Physics A, vol. 540, no. 3,
pp. 501 – 519, 1992.
[60] W. Jaus and G. Rasche, “Nuclear-structure dependence of o(α) corrections to fermi
decays and the value of the kobayashi-maskawa matrix element Vud ,” Phys. Rev. D,
vol. 41, pp. 166–176, Jan 1990.
[61] J. C. Hardy and I. S. Towner, “Superallowed 0+ → 0+ nuclear β decays: A critical
survey with tests of the conserved vector current hypothesis and the standard model,”
Phys. Rev. C, vol. 71, p. 055501, May 2005.
[62] J. Suhonen and O. Civitarese, “Weak-interaction and nuclear-structure aspects of nuclear double beta decay,” Physics Reports, vol. 300, no. 3, pp. 123 – 214, 1998.
[63] F. T. Avignone, S. R. Elliott, and J. Engel, “Double beta decay, majorana neutrinos,
and neutrino mass,” Rev. Mod. Phys., vol. 80, pp. 481–516, Apr 2008.
[64] S. R. Elliott, A. A. Hahn, and M. K. Moe, “Direct evidence for two-neutrino doublebeta decay in 82 Se,” Phys. Rev. Lett., vol. 59, pp. 2020–2023, Nov 1987.
[65] H. S. Miley, F. T. Avignone, R. L. Brodzinski, J. I. Collar, and J. H. Reeves, “Suggestive evidence for the two-neutrino double-β decay of 76 Ge,” Phys. Rev. Lett., vol. 65,
pp. 3092–3095, Dec 1990.
[66] D. L. Lincoln, J. D. Holt, G. Bollen, M. Brodeur, S. Bustabad, J. Engel, S. J. Novario, M. Redshaw, R. Ringle, and S. Schwarz, “First direct double-β decay q-value
measurement of 82 Se in support of understanding the nature of the neutrino,” Phys.
Rev. Lett., vol. 110, p. 012501, Jan 2013.

152

[67] K. Gulyuz, J. Ariche, G. Bollen, S. Bustabad, M. Eibach, C. Izzo, S. J. Novario,
M. Redshaw, R. Ringle, R. Sandler, S. Schwarz, and A. A. Valverde, “Determination
of the direct double-β-decay q value of 96 Zr and atomic masses of 90−92,94,96 Zr and
92,94−98,100 Mo,” Phys. Rev. C, vol. 91, p. 055501, May 2015.
[68] M. Redshaw, G. Bollen, M. Brodeur, S. Bustabad, D. L. Lincoln, S. J. Novario,
R. Ringle, and S. Schwarz, “Atomic mass and double-β-decay q value of 48 Ca,” Phys.
Rev. C, vol. 86, p. 041306, Oct 2012.
[69] S. Bustabad, G. Bollen, M. Brodeur, D. L. Lincoln, S. J. Novario, M. Redshaw,
R. Ringle, S. Schwarz, and A. A. Valverde, “First direct determination of the 48 Ca
double-β decay q value,” Phys. Rev. C, vol. 88, p. 022501, Aug 2013.
[70] K. Hirata, T. Kajita, M. Koshiba, M. Nakahata, Y. Oyama, N. Sato, A. Suzuki,
M. Takita, Y. Totsuka, T. Kifune, T. Suda, K. Takahashi, T. Tanimori, K. Miyano,
M. Yamada, E. W. Beier, L. R. Feldscher, S. B. Kim, A. K. Mann, F. M. Newcomer,
R. Van, W. Zhang, and B. G. Cortez, “Observation of a neutrino burst from the
supernova sn1987a,” Phys. Rev. Lett., vol. 58, pp. 1490–1493, Apr 1987.
[71] R. M. Bionta, G. Blewitt, C. B. Bratton, D. Casper, A. Ciocio, R. Claus, B. Cortez,
M. Crouch, S. T. Dye, S. Errede, G. W. Foster, W. Gajewski, K. S. Ganezer, M. Goldhaber, T. J. Haines, T. W. Jones, D. Kielczewska, W. R. Kropp, J. G. Learned, J. M.
LoSecco, J. Matthews, R. Miller, M. S. Mudan, H. S. Park, L. R. Price, F. Reines,
J. Schultz, S. Seidel, E. Shumard, D. Sinclair, H. W. Sobel, J. L. Stone, L. R. Sulak,
R. Svoboda, G. Thornton, J. C. van der Velde, and C. Wuest, “Observation of a neutrino burst in coincidence with supernova 1987a in the large magellanic cloud,” Phys.
Rev. Lett., vol. 58, pp. 1494–1496, Apr 1987.
[72] C. Hargrove, I. Batkin, M. Sundaresan, and J. Dubeau, “A lead astronomical neutrino
detector: Land,” Astroparticle Physics, vol. 5, no. 2, pp. 183 – 196, 1996.
[73] D. B. Cline, G. M. Fuller, W. P. Hong, B. Meyer, and J. Wilson, “Prospects for detection of a cosmologically significant neutrino mass from a galactic supernova neutrino
burst using a neutral-current-based detector,” Phys. Rev. D, vol. 50, pp. 720–729, Jul
1994.
[74] G. Ewan, “The sudbury neutrino observatory,” Nuclear Instruments and Methods in
Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated
Equipment, vol. 314, no. 2, pp. 373 – 379, 1992.
[75] K. Langanke, P. Vogel, and E. Kolbe, “Signal for supernova νµ and ντ neutrinos in
water Čerenkov detectors,” Phys. Rev. Lett., vol. 76, pp. 2629–2632, Apr 1996.
[76] E. Kolbe, K. Langanke, G. Martnez-Pinedo, and P. Vogel, “Neutrinonucleus reactions
and nuclear structure,” Journal of Physics G: Nuclear and Particle Physics, vol. 29,
no. 11, p. 2569, 2003.

153

[77] O. Benhar, N. Farina, H. Nakamura, M. Sakuda, and R. Seki, “Electron- and neutrinonucleus scattering in the impulse approximation regime,” Phys. Rev. D, vol. 72,
p. 053005, Sep 2005.
[78] E. Fermi, “Versuch einer Theorie der β-Strahlen. I,” Zeitschrift fur Physik, vol. 88,
pp. 161–177, Mar. 1934.
[79] F. L. Wilson, “Fermi’s theory of beta decay,” American Journal of Physics, vol. 36,
no. 12, pp. 1150–1160, 1968.
[80] B. H. Wildenthal, M. S. Curtin, and B. A. Brown, “Predicted features of the beta
decay of neutron-rich sd-shell nuclei,” Phys. Rev. C, vol. 28, pp. 1343–1366, Sep 1983.
[81] B. A. Brown and B. H. Wildenthal, “Corrections to the free-nucleon values of the singleparticle matrix elements of the m1 and gamow-teller operators, from a comparison of
shell-model predictions with sd-shell data,” Phys. Rev. C, vol. 28, pp. 2397–2413, Dec
1983.
[82] E. K. Warburton and B. A. Brown, “Effective interactions for the 0p1s0d nuclear
shell-model space,” Phys. Rev. C, vol. 46, pp. 923–944, Sep 1992.
[83] W. E. Ormand and B. A. Brown, “Isospin-mixing corrections for fp-shell fermi transitions,” Phys. Rev. C, vol. 52, pp. 2455–2460, Nov 1995.
[84] I. S. Towner and F. C. Khanna, “Quenching of allowed gamow-teller β transitions in
mirror nuclei,” Phys. Rev. Lett., vol. 42, pp. 51–54, Jan 1979.
[85] A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, S. Bacca, B. Carlsson, and D. Gazit, “Effects of three-nucleon forces and two-body currents on gamowteller strengths,” Phys. Rev. Lett., vol. 113, p. 262504, Dec 2014.
[86] S. K. Bogner, H. Hergert, J. D. Holt, A. Schwenk, S. Binder, A. Calci, J. Langhammer,
and R. Roth, “Nonperturbative shell-model interactions from the in-medium similarity
renormalization group,” Phys. Rev. Lett., vol. 113, p. 142501, Oct 2014.
[87] J. Suhonen and O. Civitarese, “Probing the quenching of ga by single and double beta
decays,” Physics Letters B, vol. 725, no. 1, pp. 153 – 157, 2013.
[88] J. Engel, “Nuclear matrix elements for double-b decay,” XXVI International Conference on Neutrino Physics and Astrophysics (Neutrino 2014), vol. 1666, 2015.
[89] J. C. Slater, “The theory of complex spectra,” Phys. Rev., vol. 34, pp. 1293–1322, Nov
1929.
[90] A. Fetter and J. Walecka, Quantum Theory of Many-particle Systems. Dover Books
on Physics, Dover Publications, 2003.
[91] G. C. Wick, “The evaluation of the collision matrix,” Phys. Rev., vol. 80, pp. 268–272,
Oct 1950.

154

[92] D. R. Hartree, “The wave mechanics of an atom with a non-coulomb central field.
part i. theory and methods,” Mathematical Proceedings of the Cambridge Philosophical
Society, vol. 24, no. 1, p. 89110, 1928.
[93] V.
Fock,
“Näherungsmethode
zur
lösung
des
quantenmechanischen
mehrkörperproblems,” Zeitschrift für Physik, vol. 61, pp. 126–148, Jan 1930.
[94] D. Thouless, “Stability conditions and nuclear rotations in the hartree-fock theory,”
Nuclear Physics, vol. 21, pp. 225 – 232, 1960.
[95] L. Brillouin, “Les problèmes de perturbations et les champs self-consistents,” Le journal
de physique et le radium, vol. 3, no. 9, pp. 373–389, 1932.
[96] P. Pulay, “Convergence acceleration of iterative sequences. the case of scf iteration,”
Chem. Phys. Lett., vol. 73, no. 2, pp. 393 – 398, 1980.
[97] P. Pulay, “Improved scf convergence acceleration,” J. Comput. Chem., vol. 3, no. 4,
pp. 556–560, 1982.
[98] C. G. Broyden, “A class of methods for solving nonlinear simultaneous equations,”
Math. Comput., vol. 19, no. 92, pp. 577–593, 1965.
[99] E. U. Condon, “The theory of complex spectra,” Phys. Rev., vol. 36, pp. 1121–1133,
Oct 1930.
[100] C. D. Sherrill and H. F. Schaefer, “The configuration interaction method: Advances
in highly correlated approaches,” vol. 34 of Advances in Quantum Chemistry, pp. 143
– 269, Academic Press, 1999.
[101] C. Møller and M. S. Plesset, “Note on an approximation treatment for many-electron
systems,” Phys. Rev., vol. 46, pp. 618–622, Oct 1934.
[102] J. Rayleigh, The Theory of Sound. The Theory of Sound, Macmillan, 1894.
[103] E. Schrödinger, “Quantification of the eigen value problem,” Annalen der Physik,
vol. 80, no. 13, pp. 437–490, 1926.
[104] L. M. Frantz and R. L. Mills, “Many-body basis for the optical model,” Nuclear
Physics, vol. 15, pp. 16 – 32, 1960.
[105] B. H. Brandow, “Linked-cluster expansions for the nuclear many-body problem,” Rev.
Mod. Phys., vol. 39, pp. 771–828, Oct 1967.
[106] T. P. Živković, “Existence and reality of solutions of the coupled-cluster equations,”
International Journal of Quantum Chemistry, vol. 12, no. S11, pp. 413–420, 1977.
[107] P. Piecuch, S. Zarrabian, J. Paldus, and J. Čı́žek, “Coupled-cluster approaches with an
approximate account of triexcitations and the optimized-inner-projection technique.
ii. coupled-cluster results for cyclic-polyene model systems,” Phys. Rev. B, vol. 42,
pp. 3351–3379, Aug 1990.
155

[108] P. Piecuch and K. Kowalski, “In search of the relationship between multiple solutions characterizing coupled-cluster theories,” in Computational Chemistry: Reviews
of Current Trends (J. Leszczynski, ed.), vol. 5, pp. 1–104, Singapore: World Scientific,
2000.
[109] N. Moiseyev, Non-Hermitian Quantum Mechanics. Cambridge University Press, 2011.
[110] R. J. Bartlett and M. Musial, “Coupled-cluster theory in quantum chemistry,” Rev.
Mod. Phys., vol. 79, pp. 291–352, Feb 2007.
[111] G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, “Coupled-cluster computations of atomic nuclei,” Rep. Prog. in Phys., vol. 77, p. 096302, Sep 2014.
[112] P. Szakács and P. R. Surján, “Stability conditions for the coupled cluster equations,”
International Journal of Quantum Chemistry, vol. 108, no. 12, pp. 2043–2052, 2008.
[113] “An updated set of basic linear algebra subprograms (blas),” ACM Trans. Math.
Softw., vol. 28, p. 135, 2002.
[114] E. K. U. Gross, E. Runge, and O. Heinonen, Many-particle theory / E.K.U. Gross, E.
Runge, O. Heinonen. A. Hilger Bristol ; Philadelphia, 1991.
[115] N. Michel, W. Nazarewicz, M. Ploszajczak, and J. Rotureau, “Antibound states and
halo formation in the gamow shell model,” Phys. Rev. C, vol. 74, p. 054305, Nov 2006.
[116] G. Hagen, D. Dean, M. Hjorth-Jensen, and T. Papenbrock, “Complex coupled-cluster
approach to an ab-initio description of open quantum systems,” Physics Letters B,
vol. 656, no. 4, pp. 169 – 173, 2007.
[117] E. Wigner and J. Griffin, Group Theory and Its Application to the Quantum Mechanics
of Atomic Spectra. Pure and applied Physics, Academic Press, 1959.
[118] C. Eckart, “The application of group theory to the quantum dynamics of monatomic
systems,” Rev. Mod. Phys., vol. 2, pp. 305–380, Jul 1930.
[119] T. Hatsuda and T. Kunihiro, “Qcd phenomenology based on a chiral effective lagrangian,” Physics Reports, vol. 247, no. 5, pp. 221 – 367, 1994.
[120] G. P. Lepage and S. J. Brodsky, “Exclusive processes in perturbative quantum chromodynamics,” Phys. Rev. D, vol. 22, pp. 2157–2198, Nov 1980.
[121] S. R. Beane, E. Chang, W. Detmold, H. W. Lin, T. C. Luu, K. Orginos, A. Parreño,
M. J. Savage, A. Torok, and A. Walker-Loud, “Deuteron and exotic two-body bound
states from lattice qcd,” Phys. Rev. D, vol. 85, p. 054511, Mar 2012.
[122] J. Beringer, J. F. Arguin, R. M. Barnett, K. Copic, O. Dahl, D. E. Groom, C. J.
Lin, J. Lys, H. Murayama, C. G. Wohl, W. M. Yao, P. A. Zyla, C. Amsler, M. Antonelli, D. M. Asner, H. Baer, H. R. Band, T. Basaglia, C. W. Bauer, J. J. Beatty,
V. I. Belousov, E. Bergren, G. Bernardi, W. Bertl, S. Bethke, H. Bichsel, O. Biebel,

156

E. Blucher, S. Blusk, G. Brooijmans, O. Buchmueller, R. N. Cahn, M. Carena, A. Ceccucci, D. Chakraborty, M. C. Chen, R. S. Chivukula, G. Cowan, G. D’Ambrosio,
T. Damour, D. de Florian, A. de Gouvêa, T. DeGrand, P. de Jong, G. Dissertori,
B. Dobrescu, M. Doser, M. Drees, D. A. Edwards, S. Eidelman, J. Erler, V. V. Ezhela,
W. Fetscher, B. D. Fields, B. Foster, T. K. Gaisser, L. Garren, H. J. Gerber, G. Gerbier, T. Gherghetta, S. Golwala, M. Goodman, C. Grab, A. V. Gritsan, J. F. Grivaz, M. Grünewald, A. Gurtu, T. Gutsche, H. E. Haber, K. Hagiwara, C. Hagmann,
C. Hanhart, S. Hashimoto, K. G. Hayes, M. Heffner, B. Heltsley, J. J. HernándezRey, K. Hikasa, A. Höcker, J. Holder, A. Holtkamp, J. Huston, J. D. Jackson, K. F.
Johnson, T. Junk, D. Karlen, D. Kirkby, S. R. Klein, E. Klempt, R. V. Kowalewski,
F. Krauss, M. Kreps, B. Krusche, Y. V. Kuyanov, Y. Kwon, O. Lahav, J. Laiho, P. Langacker, A. Liddle, Z. Ligeti, T. M. Liss, L. Littenberg, K. S. Lugovsky, S. B. Lugovsky,
T. Mannel, A. V. Manohar, W. J. Marciano, A. D. Martin, A. Masoni, J. Matthews,
D. Milstead, R. Miquel, K. Mönig, F. Moortgat, K. Nakamura, M. Narain, P. Nason,
S. Navas, M. Neubert, P. Nevski, Y. Nir, K. A. Olive, L. Pape, J. Parsons, C. Patrignani, J. A. Peacock, S. T. Petcov, A. Piepke, A. Pomarol, G. Punzi, A. Quadt, S. Raby,
G. Raffelt, B. N. Ratcliff, P. Richardson, S. Roesler, S. Rolli, A. Romaniouk, L. J.
Rosenberg, J. L. Rosner, C. T. Sachrajda, Y. Sakai, G. P. Salam, S. Sarkar, F. Sauli,
O. Schneider, K. Scholberg, D. Scott, W. G. Seligman, M. H. Shaevitz, S. R. Sharpe,
M. Silari, T. Sjöstrand, P. Skands, J. G. Smith, G. F. Smoot, S. Spanier, H. Spieler,
A. Stahl, T. Stanev, S. L. Stone, T. Sumiyoshi, M. J. Syphers, F. Takahashi, M. Tanabashi, J. Terning, M. Titov, N. P. Tkachenko, N. A. Törnqvist, D. Tovey, G. Valencia,
K. van Bibber, G. Venanzoni, M. G. Vincter, P. Vogel, A. Vogt, W. Walkowiak, C. W.
Walter, D. R. Ward, T. Watari, G. Weiglein, E. J. Weinberg, L. R. Wiencke, L. Wolfenstein, J. Womersley, C. L. Woody, R. L. Workman, A. Yamamoto, G. P. Zeller, O. V.
Zenin, J. Zhang, R. Y. Zhu, G. Harper, V. S. Lugovsky, and P. Schaffner, “Review of
particle physics,” Phys. Rev. D, vol. 86, p. 010001, Jul 2012.
[123] D. R. Entem and R. Machleidt, “Accurate charge-dependent nucleon-nucleon potential
at fourth order of chiral perturbation theory,” Phys. Rev. C, vol. 68, p. 041001, Oct
2003.
[124] P. Navrátil, “Local three-nucleon interaction from chiral effective field theory,” FewBody Systems, vol. 41, pp. 117–140, Dec 2007.
[125] R. Roth, S. Binder, K. Vobig, A. Calci, J. Langhammer, and P. Navrátil, “Mediummass nuclei with normal-ordered chiral N N + 3N interactions,” Phys. Rev. Lett.,
vol. 109, p. 052501, Jul 2012.
[126] S. Levit and D. Orgad, “Statistics of hartree-fock levels in small disordered systems,”
Physical review. B, Condensed matter, vol. 60, no. 8, pp. 5549–5560, 1999.
[127] H. J. Lipkin, “Center-of-mass motion in the nuclear shell model,” Phys. Rev., vol. 110,
pp. 1395–1397, Jun 1958.
[128] D. Gloeckner and R. Lawson, “Spurious center-of-mass motion,” Physics Letters B,
vol. 53, no. 4, pp. 313 – 318, 1974.
157

[129] R. F. Bishop, M. F. Flynn, M. C. Boscá, E. Buenda, and R. Guardiola, “Translationally
invariant coupled cluster theory for simple finite systems,” Phys. Rev. C, vol. 42,
pp. 1341–1360, Oct 1990.
[130] A. Nogga, H. Kamada, W. Glöckle, and B. R. Barrett, “The α particle based on
modern nuclear forces,” Phys. Rev. C, vol. 65, p. 054003, May 2002.
[131] C. M. Vincent, “Optimal separation of center-of-mass motion,” Phys. Rev. C, vol. 8,
pp. 929–937, Sep 1973.
[132] G. Hagen, T. Papenbrock, and D. J. Dean, “Solution of the center-of-mass problem in
nuclear structure calculations,” Phys. Rev. Lett., vol. 103, p. 062503, Aug 2009.
[133] G. R. Jansen, “Spherical coupled-cluster theory for open-shell nuclei,” Phys. Rev. C,
vol. 88, p. 024305, Aug 2013.
[134] G. H. Diercksen, B. O. Roos, and A. J. Sadlej, “Legitimate calculation of first-order
molecular properties in the case of limited ci functions. dipole moments,” Chemical
Physics, vol. 59, no. 1, pp. 29 – 39, 1981.
[135] M. Ernzerhof, C. Marian, and S. Peyerimhoff, “Energy derivative versus expectation
value approach: the dipole moment of co,” Chemical Physics Letters, vol. 204, no. 1,
pp. 59 – 64, 1993.
[136] D. J. Rowe, “Equations-of-motion method and the extended shell model,” Rev. Mod.
Phys., vol. 40, pp. 153–166, Jan 1968.
[137] C. Lanczos, “An iteration method for the solution of the eigenvalue problem of linear
differential and integral operators,” Journal of research of the National Bureau of
Standards (1934), vol. 45, pp. 255–282, 1950 10.
[138] W. E. Arnoldi, “The principle of minimized iterations in the solution of the matrix
eigenvalue problem,” Quarterly of applied mathematics, vol. 9, no. 1, pp. 17–29, 1951.
[139] R. Lehoucq, D. Sorensen, and C. Yang, ARPACK Users’ Guide. Society for Industrial
and Applied Mathematics, 1998.
[140] S. M. Reimann and M. Manninen, “Electronic structure of quantum dots,” Rev. Mod.
Phys., vol. 74, p. 1238, 2002.
[141] H.-A. Engel, V. N. Golovach, D. Loss, L. M. K. Vandersypen, J. M. Elzerman, R. Hanson, and L. P. Kouwenhoven, “Measurement Efficiency and n-Shot Readout of Spin
Qubits,” Phys. Rev. Lett., vol. 93, p. 106804, 2004.
[142] J. Birman, R. Nazmitdinov, and V. Yukalov, “Effects of symmetry breaking in finite
quantum systems,” Phys. Rep., vol. 526, no. 1, pp. 1 – 91, 2013. Effects of symmetry
breaking in finite quantum systems.

158

[143] S. Tarucha, D. G. Austing, T. Honda, R. J. van der Hage, and L. P. Kouwenhoven,
“Shell Filling and Spin Effects in a Few Electron Quantum Dot,” Phys. Rev. Lett.,
vol. 77, p. 3613, 1996.
[144] M. P. Lohne, “Coupled-cluster studies of quantum dots,” Master’s thesis, University
of Oslo, 2010.
[145] “NIST Digital Library of Mathematical Functions.” http://dlmf.nist.gov/, Release
1.0.14 of 2016-12-21, 2016. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I.
Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, eds.
[146] I. Lindgren, “The rayleigh-schrodinger perturbation and the linked-diagram theorem
for a multi-configurational model space,” J. Phys. Pt. B Atom. M. P., vol. 7, no. 18,
p. 2441, 1974.
[147] J. Høgberget, “Quantum Monte-Carlo Studies of Generalized Many-body Systems,”
Master’s thesis, University of Oslo, 2013.
[148] V. K. B. Olsen, “Full Configuration Interaction Simulation of Quantum Dots,” Master’s
thesis, University of Oslo, 2013.
[149] N. M. Parzuchowski, T. D. Morris, and S. K. Bogner, “Ab initio excited states from
the in-medium similarity renormalization group,” Phys. Rev. C, vol. 95, p. 044304,
Apr 2017.
[150] K. Emrich, “An extension of the coupled cluster formalism to excited states (i),”
Nuclear Physics A, vol. 351, no. 3, pp. 379 – 396, 1981.
[151] K. Emrich, “An extension of the coupled cluster formalism to excited states: (ii).
approximations and tests,” Nuclear Physics A, vol. 351, no. 3, pp. 397 – 438, 1981.
[152] K. Emrich and J. Zabolitzky, “Negative parity states in 16 O from coupled-cluster
equations,” Nuclear Physics A, vol. 351, no. 3, pp. 439 – 455, 1981.
[153] J. F. Stanton and R. J. Bartlett, “The equation of motion coupledcluster method. a
systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties,” The Journal of Chemical Physics, vol. 98, no. 9,
pp. 7029–7039, 1993.
[154] G. Hagen, M. Hjorth-Jensen, G. R. Jansen, R. Machleidt, and T. Papenbrock, “Continuum effects and three-nucleon forces in neutron-rich oxygen isotopes,” Phys. Rev.
Lett., vol. 108, p. 242501, Jun 2012.
[155] A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, B. D. Carlsson,
C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, “Accurate nuclear radii
and binding energies from a chiral interaction,” Phys. Rev. C, vol. 91, p. 051301, May
2015.

159

[156] F. Halzen and A. D. Martin, Quarks and Leptons: An Introductory Course in Modern
Particle Physics. 1984.
[157] J. C. Hardy and I. S. Towner, “Superallowed 0+ → 0+ nuclear β decays: A critical
survey with tests of the conserved vector current hypothesis and the standard model,”
Phys. Rev. C, vol. 71, p. 055501, May 2005.
[158] I. S. Towner and J. C. Hardy, “Improved calculation of the isospin-symmetry-breaking
corrections to superallowed fermi β decay,” Phys. Rev. C, vol. 77, p. 025501, Feb 2008.
[159] K. Kubodera, J. Delorme, and M. Rho, “Axial currents in nuclei,” Phys. Rev. Lett.,
vol. 40, pp. 755–758, Mar 1978.
[160] E. Oset and M. Rho, “Axial currents in nuclei: The gamow-teller matrix element,”
Phys. Rev. Lett., vol. 42, pp. 47–50, Jan 1979.
[161] D. Wilkinson, “Analysis of neutron beta decay,” Nuclear Physics A, vol. 377, no. 2,
pp. 474 – 504, 1982.
[162] J. Suhonen, From Nucleons to Nucleus. Theoretical and Mathematical Physics, Berlin,
Germany: Springer, 2007.

160