COMPUTATIONAL DEVELOPMENTS FOR AB INITIO MANY-BODY THEORY
By
Justin Gage Lietz
A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
Physics - Doctor of Philosophy
Computational Math, Science and Engineering - Dual Major
2019
COMPUTATIONAL DEVELOPMENTS FOR AB INITIO MANY-BODY
ABSTRACT
THEORY
By
Justin Gage Lietz
Quantum many-body physics is the body of knowledge which studies systems of many
interacting particles and the mathematical framework for calculating properties of these
systems. Methods in many-body physics which use a ﬁrst principles approach to solving
the many-body Schr¨odinger equation are referred to as ab initio methods, and provide ap-
proximate solutions which are systematically improvable. Coupled cluster theory is an ab
initio quantum many-body method which has been shown to provide accurate calculations of
ground state energies for a wide range of systems in quantum chemistry and nuclear physics.
Calculations of physical properties using ab initio many-body methods can be computa-
tionally expensive, requiring the development of eﬃcient data structures, algorithms and
techniques in high-performance computing to achieve numerical accuracy.
Many physical systems of interest are diﬃcult or impossible to measure experimentally,
and so are reliant on predictive and accurate calculations from many-body theory. Neutron
stars in particular are diﬃcult to collect observational data for, but simulations of inﬁnite
nuclear matter can provide key insights to the internal structure of these astronomical ob-
jects. The main focus of this thesis is the development of a large and versatile coupled cluster
program which implements a sparse tensor storage scheme and eﬃcient tensor contraction
algorithms. A distributed memory data structure for these large, sparse tensors is used so
that the code can run in a high-performance computing setting, and can thus handle the
computational challenges of inﬁnite nuclear matter calculations using large basis sets. By
validating these data structures and algorithms in the context of coupled cluster theory and
inﬁnite nuclear matter, they can be applied to a wide range of many-body methods and
physical systems.
Dedicated to my loving parents.
iv
ACKNOWLEDGMENTS
This dissertation is the culmination of six years of graduate school, all of which has been
under the guidance of my PhD advisor, Morten Hjorth-Jensen.
I am incredibly lucky to
have been on this journey with Morten, as he is not only an excellent teacher, scientist, and
mentor, but he is also a wonderful person who is overﬂowing with kindness. It is hard to
imagine what this dissertation would look like if I had never met Morten, as his inﬂuence is
in every one of its pages. My thoughts surrounding physics, research, and philosophy have
been profoundly shaped by Morten and that extends much beyond just this work. Thank
you so much Morten, while this dissertation marks the end of my time as your student, I
know it also marks the beginning of a lifelong friendship.
Next, I would like to thank Scott Bogner and Heiko Hergert who also have had a huge
impact on this dissertation by creating an intellectually stimulating and friendly environment
for the nuclear many-body group. It has been a pleasure to share this journey with you and
with the other grad students in the group: Sam Novario, Nathan Parzuchowski, Titus Morris
and Fei Yuan. The times I got to travel to summer schools with all of you and explore new
cities are some of my favorite memories of grad school.
Outside of those who directly contributed to this dissertion, I also want to thank those
who have made a positive impact on my life while I was working on it. First I want to thank
my girlfriend, Rachel, for all of your love and support. You helped make my time in grad
school the best years of my life, and I can’t wait for our next adventure together. I am also
grateful for the time shared with my other grad student friends, and in particular the great
times I got to spend with Dennis, Terri, Michael, and Brent in the last few months. It was
tough to push through at the end of grad school, but getting to spend time with you all
v
helped immensely.
Next I would like to thank my brother, Jordan who has been my friend and role model
for my entire life and continues to be to this day, and my sister-in-law Sarah who can bring
a smile to anyone in my family. You guys are the best, and I can’t wait for many more visits
to see you, Annabelle, and Ollie. Lastly, and certainly most importantly, I have to give
thanks to my parents, Kim and Steve, to whom this dissertation is dedicated to. They have
provided unconditional love and selﬂess support so that I could pursue what I love doing,
and I am incredibly grateful for everything they have done for me.
Thank you so much mom and dad.
vi
TABLE OF CONTENTS
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ix
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
x
KEY TO ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii
Chapter 1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Nuclear Theory as a Window into the the Stars
. . . . . . . . . . . . . . . .
1.2 Ab initio methods in Nuclear Theory . . . . . . . . . . . . . . . . . . . . . .
1.3 Quantum Many-Body Methods
. . . . . . . . . . . . . . . . . . . . . . . . .
1.4
Inﬁnite Matter Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Computational Challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.6 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Chapter 2 Quantum Many-Body Physics . . . . . . . . . . . . . . . . . . . .
2.1 Bra-ket Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Many-fermion wave functions and spaces . . . . . . . . . . . . . . . . . . . .
2.3 Occupation Number Representation . . . . . . . . . . . . . . . . . . . . . . .
2.4 Creation and annihilation operators . . . . . . . . . . . . . . . . . . . . . . .
2.5 Number Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 Anti-commutation relations
. . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7 Operators in Second Quantized Form . . . . . . . . . . . . . . . . . . . . . .
2.8 Wick’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.9 Generalized Wick’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.10 Slater-Condon Rules
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.11 The Fermi Vacuum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.12 Conﬁguration Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.13 Hartree-Fock Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.14 Many-Body Perturbation Theory . . . . . . . . . . . . . . . . . . . . . . . .
2.15 In-Medium Similarity Renormalization Group . . . . . . . . . . . . . . . . .
2.16 The Magnus Formulation of IM-SRG . . . . . . . . . . . . . . . . . . . . . .
Chapter 3 Physical Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1 Pairing Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Single-Particle Basis for Inﬁnite Fermionic Matter . . . . . . . . . . . . . . .
3.3 Two-Nucleon Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Homogeneous Electron Gas
. . . . . . . . . . . . . . . . . . . . . . . . . . .
Chapter 4 Coupled Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1 Prologue to Coupled Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1
3
4
5
6
8
9
11
14
21
23
26
28
29
33
38
38
44
50
53
58
65
72
84
85
89
93
97
99
99
vii
4.2 Coupled Cluster Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.3 Coupled Cluster Diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
4.4 Diagrammatic Derivation of the Coupled Cluster Equations . . . . . . . . . . 117
4.5 Computational Scaling of Coupled Cluster Theory . . . . . . . . . . . . . . . 127
Inﬁnite Neutron Matter
Chapter 5 Computational Methodology . . . . . . . . . . . . . . . . . . . . . 132
5.1 Code Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
5.1.1 Pairing Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
. . . . . . . . . . . . . . . . . . . . . . . . . 135
5.1.2
5.2 Taming the Two-Body Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
5.3 Performance Testing Matrix-Matrix Multiplication . . . . . . . . . . . . . . . 141
5.4 Tensor Contractions as Matrix Multiplication . . . . . . . . . . . . . . . . . 145
5.5 Parallel Computing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
5.6 Distributed Memory Parallelization . . . . . . . . . . . . . . . . . . . . . . . 156
5.7 Final Parallel Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
Chapter 6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
6.1 Neutron Matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
6.2 Homogeneous Electron Gas
. . . . . . . . . . . . . . . . . . . . . . . . . . . 171
6.3 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
Chapter 7 Conclusions and Perspectives . . . . . . . . . . . . . . . . . . . . . 184
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
viii
LIST OF TABLES
Table 2.1:
Single-Particle Index Conventions
. . . . . . . . . . . . . . . . . . .
45
Table 3.1:
Single-particle states and their quantum numbers and their energies
from Eq. (3.3). The degeneracy for every quantum number p is equal
. . . . . . . . . . . . . .
to two due to the two possible spin values.
Table 3.2:
y + n2
Total number of particle ﬁlling N↑↓ for various n2
z values
for one spin-1/2 fermion species. Borrowing from nuclear shell-model
terminology, ﬁlled shells correspond to all single-particle states for
one n2
z value being occupied. For matter with both protons
and neutrons, the ﬁlling degree increased with a factor of 2. . . . . .
x + n2
x + n2
y + n2
86
92
Table 3.3:
Parameters used to deﬁne the Minnesota interaction model [53].
. .
97
Table 4.1:
Sign table for the four terms of eq. (4.37). . . . . . . . . . . . . . . . 115
Table 5.1:
Coupled cluster and MBPT2 results for the simple pairing model with
eight single-particle levels and four spin-1/2 fermions for diﬀerent
values of the interaction strength g.
. . . . . . . . . . . . . . . . . . 134
Table 5.2:
Table 5.3:
CCD and MBPT2 results for inﬁnite neutron matter with N = 66
neutrons and a maximum number of single-particle states constrained
by Nmax = 36 (36 plane wave energy shells).
. . . . . . . . . . . . . 136
All possible particle-hole sectors are listed in the left column. In the
right column are 6 particle hole sectors which contain all of the infor-
mation of the whole matrix, plus how the other 10 can be equivalently
expressed.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
Table 5.4:
A straight forward scheme to organize the two-body basis in columns. 139
ix
LIST OF FIGURES
Figure 2.1:
SRG with direct integration and with the Magnus expansion.
. . . .
81
Figure 2.2: Magnus SRG with the exact unitary transformation and with a BCH
expansion truncated after a ﬁxed number of terms . . . . . . . . . .
Figure 2.3: Magnus SRG with the exact unitary transformation and with a BCH
. . . . . . . . . .
expansion truncated after a ﬁxed tolerance is met
Figure 3.1: Correlation energy for the pairing model with exact diagonalization,
MBPT2 and perturbation theory to third order MBPT3 for a range of
interaction values. A canonical Hartree-Fock basis has been employed
in all MBPT calculations. . . . . . . . . . . . . . . . . . . . . . . . .
82
83
88
Figure 5.1: 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.
. . . . . . . . . . . . . . 135
Figure 5.2: The pp-pp sector of a two-body interaction matrix for a simple neu-
tron matter system with 40 single-particle states above the Fermi
level.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
Figure 5.3:
Implementation of the same mathematics can have very diﬀerent run
times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
Figure 5.4: MBPT2 contribution to the correlation for pure neutron matter with
N = 14 neutrons and periodic boundary conditions. Up to approx-
imately 1600 single-particle states have been included in the sums
over intermediate states in Eqs. (5.7) and (5.8) . . . . . . . . . . . . 148
Figure 5.5: A cartoon of how the interaction matrix might be split into work
loads for diﬀerent threads of execution for the naive storage and the
block storage schemes.
. . . . . . . . . . . . . . . . . . . . . . . . . 155
Figure 5.6: Performance of MBPT2 calculations with increasing number of MPI
ranks. The speed of the calculation is measured in s−1, the black
data are the inverse time required to ﬁnish the calculation on the
fastest rank, and the red data are the speeds of the slowest rank.
. . 157
x
Figure 5.7: The blocks are distributed to the ranks to try and keep the num-
ber of non-zero matrix elements equal among ranks. The histogram
shows that the time it takes to load these blocks is dominated by an
enormous amount of small blocks, which is ideal for load balancing.
159
Figure 5.8: The t-amplitudes are permuted as needed for the tensor contractions
which are not aligned as a matrix-matrix product. In the histogram,
the larger blocks have begun to take more of the total time relative
to the previous step. . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
Figure 5.9: The largest tensor contractions are now performed, which are already
aligned as matrix-matrix products. The O(N 3) scaling of the matrix-
matrix product causes the larger matrices to contribute signiﬁcantly
to the total processing time in this step. . . . . . . . . . . . . . . . . 161
Figure 5.10: The t-amplitudes are summed together and the correlation energy
is calculated. If the energy has not converged to the set tolerance,
another iteration of the CC equations are performed, using these new
t-amplitudes in step 2.
. . . . . . . . . . . . . . . . . . . . . . . . . 161
Figure 6.1: Two diﬀerent energy per particle plots at low densities of neutron
matter with the Minnesota potential [53] computed in the CCD ap-
proximation with 54 neutrons and an Nmax = 100 truncation (100
plane-wave energy shells), corresponding to 10754 single particle states.165
Figure 6.2: Energy per particle of pure neutron matter computed in the CCD
approximation with the Minnesota interaction model [53]. . . . . . . 166
Figure 6.3: The relative error shows how much the CCD correlation energy is
changing between subsequent calculations at diﬀerent model spaces
sizes ranging from Nmax = 10 to 70 for neutron matter with the Min-
nesota potential at density 0.2 fm−3.
. . . . . . . . . . . . . . . . . 167
Figure 6.4: Finite-size eﬀects in diﬀerent energies of pure neutron matter com-
puted with the Minnesota interaction model [53] as a function of the
number of particles for both periodic boundary conditions (PBC) and
twist-averaged boundary conditions (TABC5).
. . . . . . . . . . . . 169
Figure 6.5: Energy per particle for pure neutron matter with the Minnesota po-
tential [53]. Here the calculations have been performed with IM-
SRG(2), CCD, CIMC [77], and the ADC(3) Self-Consistent Green’s
Function scheme [77].
. . . . . . . . . . . . . . . . . . . . . . . . . . 170
Figure 6.6: The CCD energy per particle for the homogeneous electron gas for a
range of Wigner-Seitz Radii with A = 14 electrons . . . . . . . . . . 172
xi
Figure 6.7: Contributions to the energy from purely CCD many-body correlations.173
Figure 6.8: Fractional contribution to the energy from the Hartree-Fock reference
state.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
Figure 6.9: The relative error shows how much the CCD correlation energy is
changing between model spaces sizes ranging from Nmax = 10 to 60
for the electron gas at rs = 0.5. . . . . . . . . . . . . . . . . . . . . . 174
Figure 6.10: The relative error shows how much the CCD correlation energy is
changing between model spaces sizes ranging from Nmax = 10 to 60
for the electron gas at rs = 0.1. . . . . . . . . . . . . . . . . . . . . . 175
Figure 6.11: The relative error of the CCD correlation energy is changing between
model spaces sizes ranging from Nmax = 100 to 200 for the neutron
matter with the Minnesota potential at A = 54 and ρ = 0.08. . . . . 176
Figure 6.12: The time required for the large basis set Minnesota potential calcu-
lations.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
Figure 6.13: Strong scaling of distributed memory code, dark green line shows
ideal case.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
Figure 6.14: Weak scaling of distributed memory code, the dark green line shows
the ideal case.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
Figure 6.15: Number of tensor elements required for the 3-body force in the inﬁnite
matter basis with and without block-diagonal compression.
. . . . . 181
Figure 6.16: Size of tensor in gigabytes required for the 3-body force in the inﬁnite
matter basis with and without block-diagonal compression.
. . . . . 182
Figure 6.17: On node timing tests for the tensor contraction of three-body force
diagrams. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
xii
KEY TO ABBREVIATIONS
NSCL - National Superconducting Cyclotron Laboratory (East Lansing MI, USA)
CC - Coupled Cluster
CC(S)(D)(T) - Coupled Cluster with (Singles)(Doubles)(Triples)
IM-SRG - In-Medium Similarity Renormalization Group
EoS - Equation of State
d.o.f. - Degree of Freedom
QCD - Quantum Chromodynamics
LQCD - Lattice QCD
GPU - Graphics Processing Unit
χ-EFT - Chiral Eﬀective Field Theory
FCI - Full Conﬁguration Interation
CISD - Conﬁguration Interaction Singles and Doubles
CIMC - Conﬁguration Interaction Monte Carlo
MBPT - Many-Body Perturbation Theory
BCH - Baker - Campbell - Hausdorﬀ
HEG - Homogeneous Electron Gas
BLAS - Basic Linear Algebra Subroutines
cuBLAS - CUDA BLAS
OMP - Open Multi-Processing
MPI - Message Passing Interface
ADC - Algebraic Diagrammatic Construction
xiii
Chapter 1
Introduction
1.1 Nuclear Theory as a Window into the the Stars
Stars are objects of extreme forces. Spheres of crushing gravitational forces being held up
by the violent nuclear reactions at their core. Until they aren’t. Eventually the nuclear fuel
at the core of every star depletes, giving way to gravitational contraction, and in some cases
cataclysmic collapse. Many bright burning stars will collapse onto their own cores, creating
one of the brightest events in the galaxy, a supernova. The extreme explosions of some dying
stars blast oﬀ a signiﬁcant fraction of their total material, and leave behind a spent core.
If the resulting core is large enough, the gravitational compression will crush the remaining
matter to extreme densities. In the cases that are not quite large enough for a black hole to
form, a neutron star remains at the center of a once bright burning star. Neutron stars get
their name because the extreme densities caused by the gravitational collapse have pushed
beyond the limits of electron degeneracy, collapsing the bulk of protons and electrons into
neutrons. In a region around 1-3 solar masses, the resulting neutron degeneracy pressure
together with the very close range nuclear force are enough to push back against further
gravitational collapse, forming an incredible astronomical object that is composed of a very
unique state of matter. Neutron stars tend be around one solar mass, but are only about
10km in size, as the entire macroscopic object is around nuclear density. After this bright
supernova, the resulting neutron star core is left cold and dim. This makes direct observation
1
of these fascinating objects very diﬃcult, and can only rarely be done. In the case of pulsars,
direct electromagnetic radiation can be detected, however it is only from a directed beam
out of the magnetic axis of the star. Otherwise, indirect measurements of nearby stars must
be done. Looking for the gravitational footprint on the orbits of stars, such as a bright star
in a binary system with a neutron star can determine with some accuracy the mass of the
binary neutron star partner. A direct astronomical measurement of the neutron star radius
might not be possible, as they are just too dim and too far away. Since telescopes are largely
unable to see neutron stars, it is up to theoretical physics to help paint a picture of these
extreme objects.
To ﬁgure out how neutron stars respond to gravitational compression, and how they
eventually equilibrate to some radius, the equation of state (EoS) of the neutron star must
be known. Currently, there are many proposed equations of state [1, 2], which lead to a large
spread of possible radii given a particular observed mass. This is due to the diﬃcult nature
of calculating the nuclear EoS, which in principle requires knowing the exact composition of
a neutron star, and calculating the energetic state of this quantum system. A large part of
a neutron star is thought to be pure neutron matter, although some amount of proton and
lepton matter is likely present in a state of β-stable matter. Some theories posit that the
extreme densities towards the core of the neutron star could cause the formation of hyperons
in the nuclear matter, or that the matter could be pushed into a state of pure quarks [1].
Regardless of the composition, the resulting calculation is a diﬃcult problem of quantum
many-body physics. A comprehensive treatment of this problem would involve calculating
the strong interactions of an absurd number of particles, and many attempts have been made
using a slew of diﬀerent approximations. Regardless of the exact framework, it seems that
theoretical nuclear physics, the study of some of the smallest particles, could be our best
2
tool to determine properties of these massive celestial bodies.
1.2 Ab initio methods in Nuclear Theory
Nuclear theory as a ﬁeld is currently in an exciting period of growth due to theoretical and
computational developments in ab initio methods over the last two decades. Ab initio, latin
for “from the beginning”, is a phrase used to describe work done from ﬁrst principles, which
for nuclear theory means starting from the building blocks of the atomic nucleus: protons
and neutrons. While nucleons are composite particles made up of quarks and gluons, they are
well bound (on the order of GeV) compared to the interactions between them (MeV) so they
function well as the basic degrees of freedom (d.o.f). The “hard core” nature of the nucleon-
nucleon force has made the ab initio approach to calculating nuclear properties intractable in
the past, as this leads to the coupling of high and low momentum modes, creating diﬃculties
in calculating all but the smallest nuclei. Many strategies have been used to evade this
problem, perhaps most notable are the phenomonological models used to great eﬀect with
shell model (SM) calculations. By developing nuclear interactions using nuclear data input
near regions of interest, high accuracy calculations of properties of nuclei have been made [3].
In this approach however, many nucleonic d.o.f.’s are “frozen out”, meaning that they are
ignored and some contact with the underlying physics is lost. Another approach has been
the development of similarity renormalization group (SRG) methods, which have lead to the
proper decoupling of these high and low momentum modes, leading to softer interactions that
are able to converge much faster, and make many more calculations possible [4, 5]. Potentials
generated from chiral eﬀective ﬁeld theory (χ-EFT) [6, 7], which connect the nuclear force
to the underlying symmetries of the QCD langrangian, can now be softened with SRG [4, 5]
3
leading to a class of potentials that have connections to the underlying physics, and are
tractable for calculations of nuclei. This progress has allowed properties of medium mass
nuclei to be calculated with links to fundamental forces at modest computational cost. This
has opened a whole new way of studying nuclei, as now methods that keep track of all of
the nucleonic d.o.f.’s can be used in realistic calculations.
1.3 Quantum Many-Body Methods
There are many ways of taking on realistic calculations of nuclear properties, and at
the core of many approaches is the non-relativistic many-body Schr¨odinger equation. Low
energy nuclear physicists can get away without using relativistic quantum mechanics because
the typical binding energies of nuclei is in the range of 1 MeV to 10 MeV per nucleon,
while each nucleon itself is bound together at around 1 GeV. Even with this non-relativistic
approximation, the task of calculating any many-body problem is daunting. The second
approximation that must be done is to pick a ﬁnite basis to perform the calculation in.
These basis sets are in principle inﬁnite, but we must have a ﬁnite system for our ﬁnite
computers. Within this framework, the task of a complete energy spectrum calculation of
N particles in a basis with M single-particle states, would require diagonalizing an(cid:0)M
(cid:1) sized matrix where(cid:0)M
(cid:0)M
N !(M−N )! is a factorially growing number. This means that
for all but small systems in small basis sets, this factorial growth will quickly grow beyond
(cid:1) =
N
(cid:1) by
N
N
M !
current computational power. Full conﬁguration interaction (FCI) is a method which uses
a minimal number of approximations, and computes a nearly exact energy of the system,
but at a massive computational cost [8, 9, 10, 11]. This unfavorable scaling has led to the
development of an entire industry of approximations to solving the full many-body problem.
4
Each method approaches the many-body problem in its own way, with its own series of
advantages and disadvantages. In particular, coupled cluster theory [12, 13, 14, 15, 16] has
been in use in many-body theory since the 60’s, starting with the work of Coester [17, 18]
and Kummel [19], and saw enormous success in quantum chemistry, and more recently
nuclear theory as well [20, 21]. Coupled cluster theory is centered around a way to organize
the many-body basis by grouping states in excitation “clusters” that lead to very favorable
truncations of less important terms. Truncations are almost always needed in practical
calculations, but by restoring the truncations term by term we can systematically improve
the solution and eventually restore the exact FCI answer. With this improved many-body
basis truncation, coupled cluster theory has a favorable polynomial scaling, and sacriﬁces
only a small amount of accuracy.
In quantum chemistry, CC theory has been used to
calculate molecular properties to chemical accuracy at a fraction of the computational cost
of a total FCI calculation. In nuclear physics, coupled cluster’s many-body truncation errors
are typically minimal when compared to errors in the approximation of the nuclear forces.
With this success, CC has shown to be one of the premier ab initio many-body methods in
nuclear physics.
1.4
Inﬁnite Matter Calculations
This thesis focuses on ab initio many-body calculations relevant for neutron stars. The
ultimate goal is to learn more about neutron stars while still maintaining a link to the
underlying theory of the strong force, quantum chromodynamics, while at the same time
studying the tools and approximations needed to make this possible. Simulating an entire
star at the quantum level is an impossible task. However, by studying a small periodic chunk
5
of neutron star matter, we can extract properties like the equation of state while maintaining
these important links[1].
This idea of studying a periodic box of quantum matter is not a new idea, and the so
called “inﬁnite matter” problem is an often revisited system for many-body theorists [22].
By simulating a ﬁnite number of particles in a periodic box, an approximation to an inﬁnite
space ﬁlled with these particles can be made. Inﬁnite matter is studied in quantum chemistry
[23] to simulate electrons moving in a neverending lattice of atoms, and in nuclear physics
it can be studied as a large block of neutron star matter [2].
Many-body physicists in general study the inﬁnite matter problem as a sandbox to ex-
amine their theoretical machinery. The periodic box in which inﬁnite matter is frequently
studied leads to a natural choice of basis, that of plane waves. The ﬂat periodic boundaries
that are chosen lead to a quantization of momentum modes, which give a basis of ﬁxed
momenta waves moving through the box. Despite being a theoretically convenient basis, the
plane waves can be computationally challenging, as realistic simulations of neutron matter
can need hundreds of particles and thousands of basis states to converge to a point that
resembles a true inﬁnite slab of matter.
1.5 Computational Challenges
Despite the polynomial scaling, ab initio methods like coupled cluster theory can begin
to struggle with the computational load that thousands of basis states require. It is here
that a deep inspection of the many-body tools in use must be done. Approximations in
coupled cluster theory organize the many-body basis in terms that (typically) decrease in
importance for higher excitation levels. For inﬁnite matter calculations, the ﬁrst non-zero
6
terms are that of doubles excitations, that is, exciting pairs of particles together, in what is
called coupled cluster doubles (CCD). Even in this restrictive approximation, storage of the
two-body interaction matrix scales as M 4, where the number of single-particle basis states
M can quickly grow to 103, leading to matrix sizes of 1012 elements, which is already too
large for modern computers to deal with. While these matrices are very large, by exploiting
the sparsity of this many-body basis, these matrices can be tamed to sizes that can run on a
single large computer, or a small cluster. However, studies of modern nuclear potentials have
shown that three-body forces are necessary for accurate calculations, leading to a three-body
interaction matrix that scales as M 6. Furthermore, CCD alone is often too restrictive of an
approximation, and either partial or full inclusion of triples excitations (CCDT) is necessary
for modern state-of-the-art calculations of inﬁnite nuclear matter. These complications mean
that tools from high-performance computing are necessary to meet the precision and accuracy
demands of ab initio many-body theory. First, to even store the three-body interaction
matrix which can quickly grow to hundreds of terabytes (TB) in size, distributed memory
algorithms must be used. Only by distributing these large interaction matrices across a
computational cluster can the calculation even be started. Next, the ﬂoating point operations
(FLOPs) required by a CCDT calculation with three-body forces will scale as M 9, meaning
that massively parallel algorithms are needed to get these calculations ﬁnished in a reasonable
amount of time. This massively parallel paradigm of supercomputing usually means writing
custom code for the architecture of the computer that will be used, and in the case of the
current largest computers, this means leveraging the enormous power of graphics processing
units (GPUs) to get the job done. Along with a growing demand for computational power is
a growing problem with reproducible science and portable codes. To maximize the scientiﬁc
eﬀort of computational physicists, research that is maintainable, extensible and reproducible
7
needs code that is well tested, well documented, well designed and version controlled.
1.6 Thesis Overview
This thesis will review the basics of quantum mechanics and many-body theory as it
pertains to coupled cluster calculations in Chapter 2 to establish the language and notation
needed for the subsequent chapters. To establish the context of the ﬁeld, we will go through
the relevant theory derivations, and Chapter 3 will describe the quantum systems which will
be tackled with these methods. A focus will be on coupled cluster theory, and Chapter 4 will
show how to derive CC theory using the development of diagrammatic techniques. Once this
foundation is laid, we dive into the main object of this thesis: a large and versatile computer
program which can calculate properties of a variety of quantum systems using a variety
of many-body methods. In Chapter 5 we describe the distributed memory data stuctures
and algorithms that implement these many-body methods eﬃciently in a high performance
computing setting. Lastly, the numerical results and performance testing of the program for
inﬁnite matter calculations are discussed in Chapter 6.
8
Chapter 2
Quantum Many-Body Physics
This work will focus on ab initio calculations of many interacting particles where ab
initio, meaning “from the beginning”, refers to the fact that we want these calculations to
be as fundamental as possible. By starting from the basic building blocks, we can make
accurate predictions of properties over a wide range of systems. Choosing these degrees of
freedom is a game of compromise. Degrees of freedom that are too “macroscopic” will have
limited applicability, but degrees of freedom that are too “microscopic” yield a wide range of
applicability, but coupling the very small scales up to the large scales of the system size can
become computationally impossible. The target systems of nuclear physics frequently land in
a regime where using the constituent protons and neutrons as the degrees of freedom can be
too microscopic, as calculating properties of nuclei rapidly grows too complex. Historically,
phenomenological models, like the shell model, have had success by generating eﬀective
interactions for a few valence nucleons on top of an inert closed shell core. Accurate properties
can be computed relatively quickly by “freezing out” the nucleon degrees of freedom in the
closed shell. However, the shell model interactions rely on known experimental data in the
region of interest to ﬁt the matrix elements, and so any given shell model interaction can
only function well in this limited space. Extrapolating into regions where there is little
experimental data is very challenging, as the interaction was only tuned to the region of
interest.
In the other extreme, the nucleons that make up nuclei are themselves composite particles
9
of quarks and gluons. So a truly fundamental calculation would build up the properties with
all quark and gluon degrees of freedom active. Such calculations are done in the ﬁeld of
lattice quantum chromodynamics (LQCD) [24, 25], but due to the extreme microscopic
nature of these degrees of freedom they quickly become overwhelmingly diﬃcult, leaving
only the smallest nuclear systems accessible this way .
This work operates in between these regimes, where the quarks and gluons are frozen
out, but all of the protons and neutrons remain active in the calculation. As a trade oﬀ,
the interactions between nucleons are generated from chiral eﬀective ﬁeld theory (χ-EFT)
[6, 7], which bring the symmetries from the fundamental QCD Lagrangian. While these
calculations are expensive, truncations are made to the possible conﬁgurations of nucleons
to make them feasible, and techniques in high performance computing allow relatively large
systems to be calculated. In principle, interactions ﬁt once to few nucleon data have much
larger ranges of applicability, and would accurately compute properties from small to medium
mass nuclei to inﬁnite nuclear matter.
In reality, the current state of the art predective
models need additional data like the binding energy and radii of small to medium mass
nuclei [26, 27]. In this sense, the philosophy of ab initio quantum many-body physics is to
keep all nucleons active, and while truncations are made, there is a systematically improvable
scheme for both the interaction and the many-body correlations. To provide the foundation
to make these statements concrete, this chapter will brieﬂy walk through single-particle
quantum mechanics and then survey a few quantum many-body techniques.
10
2.1 Bra-ket Notation
A concise formalism for describing quantum states is bra-ket notation, also called Dirac
notation. In bra-ket notation, a quantum state is represented as an abstract ket, |ψ(cid:105). This
notation distinguishes itself from wave mechanics where the state is written explicitly as
ψ(x), a function of R space, or matrix mechanics where the state is expanded in some
orthonormal basis, and referenced as a set of basis components (c0, c1, . . . ). The abstract
nature of the ket allows for a formalism where derivations and manipulations can be done in
an invariant way, and a choice of coordinates or a basis can be chosen at any point where it
is convenient.
There are many texts covering the formalism of bra-ket notation, so we’ll just look at a
few interesting pieces to have footing for later discussions. In bra-ket formalism, the quantum
state |ψ(cid:105) ∈ H is an element of state space, which is a abstract complex Hilbert space that is
inﬁnite dimensional and separable (i.e. can have a countable orthonormal basis).
One property of Hilbert spaces that is central to quantum mechanics is that they are
closed under linear combinations. As a consequence, superpositions of states are themselves
states in the space
|ψ(cid:105) = c1 |ψ1(cid:105) + c2 |ψ2(cid:105) .
(2.1)
Additionally, Hilbert spaces come with an inner product IP : (H,H) → C which can be
written many diﬀerent ways. However the last notation written below using the angled
brackets is where bra-ket notation gets its name
IP(|φ(cid:105) ,|ψ(cid:105)) = (cid:104)|φ(cid:105) ,|ψ(cid:105)(cid:105) ≡ (cid:104)φ|ψ(cid:105) = (cid:104)ψ|φ(cid:105)∗ .
(2.2)
11
This leads to the deﬁnition of the bra state (cid:104)φ|. For any state |φ(cid:105) ∈ H, we can associate a
linear functional (cid:104)φ| ≡ f|φ(cid:105) : H → C, where for |ψ(cid:105) ∈ H
f|φ(cid:105)(|ψ(cid:105)) = (cid:104)φ| (|ψ(cid:105)) = IP(|φ(cid:105) ,|ψ(cid:105)) = (cid:104)φ|ψ(cid:105) .
(2.3)
Only the bra-ket notation of the inner product will be used from now on, the function
argument style was just to draw attention to the fact that these bras are linear maps from
the Hilbert space to the complex numbers. The Hermitian conjugate (conjugate transpose)
is used to go from a ket state to the correspoding bra state,
|ψ(cid:105)† = (cid:104)ψ|
(cid:104)ψ|† = |ψ(cid:105) .
(2.4)
A distinct advantage to this notation is now the projection operator Pψ onto the state |ψ(cid:105)
is compactly written as
Pψ = |ψ(cid:105)(cid:104)ψ| .
(2.5)
A set of ket vectors is considered an orthonormal basis {|ψi(cid:105) ≡ |i(cid:105)}i∈N if they satisfy the
orthonormality relation
and the completeness relation
(cid:104)i|j(cid:105) = δij,∀i, j ∈ N,
(cid:88)
i∈N|i(cid:105)(cid:104)i| = 1H.
(2.6)
(2.7)
From here, we can represent the quantum state in any basis of our choosing, by applying
12
the completeness relation (2.7) above on any ket |ψ(cid:105)
∞(cid:88)
i=0
|ψ(cid:105) =
|i(cid:105)(cid:104)i|ψ(cid:105) .
(2.8)
Next, we introduce linear operators on H. An operator A maps a ket |ψ(cid:105) into a new ket
A|ψ(cid:105) = |Aψ(cid:105) ∈ H and (cid:104)ψ| A = (cid:104)A†ψ|, where A† is the adjoint of A. The adjoint is deﬁned
by the operation of the operator under the inner product
IP(|φ(cid:105) , A|ψ(cid:105)) = IP(A† |φ(cid:105) ,|ψ(cid:105)).
(2.9)
In bra-ket notation notice that (cid:104)φ|ψ(cid:105) = (cid:104)ψ|φ(cid:105)∗, so if |ψ(cid:105) is acted on by A then projected
onto (cid:104)φ| we get
(cid:104)φ| (A|ψ(cid:105)) = (cid:104)φ|Aψ(cid:105) = (cid:104)A†ψ|φ(cid:105)∗ ,
(2.10)
which deﬁnes the bra-ket notation of the adjoint. There are some technicalities regarding
the domain of A in contrast to the domain of A† which can in general cause problems, which
are detailed in[28]. If this is not a concern, as is the case in all calculations in this work,
then bra-ket notation actually provides an equivalent interpretation that the operator A acts
on the bra state, which is then projected onto by the ket state. In the case of self-adjoint
operators where A = A†, there are no worries even in bra-ket notation since
(cid:104)φ|Aψ(cid:105) = (cid:104)Aφ|ψ(cid:105) = (cid:104)φ|A|ψ(cid:105) .
(2.11)
Self-adjoint operators are of particular interest, since they correspond to physical observables
to ground the theory in reality.
13
Since this is a work in computational physics, everything must be truncated to ﬁt the
calculations on a computer. In this case, one must select a representation for the calculation
to be carried out in, and then the basis is truncated
∞(cid:88)
i=0
|ψ(cid:105) =
N(cid:88)
i=0
(cid:104)i|ψ(cid:105)|i(cid:105) ≈
(cid:104)i|ψ(cid:105)|i(cid:105) .
(2.12)
The accuracy of this approximation will be discussed in later chapters. The quantum state is
fully encoded in N coeﬃcients of orthonormal basis states and operators are fully described
by their matrix elements (cid:104)i|A|j(cid:105). In this regime, observables are represented by Hermitian
matrices where, (cid:104)i|A|j(cid:105) = (cid:104)j|A|i(cid:105)∗, which have real eigenvalues.
Let us now consider a Hamiltonian as an operator on a Hilbert space in the following way:
The eigenkets of ˆH, denoted |i(cid:105), provide an orthonormal basis for the Hilbert space. The
spectrum of allowed energy levels of the system is given by the set of eigenvalues, denoted
{εi}, solving the equation:
ˆH |i(cid:105) = εi |i(cid:105) .
(2.13)
Since ˆH is a Hermitian operator, the energy is always a real number.
2.2 Many-fermion wave functions and spaces
To calculate energies of quantum many-body systems, we must ﬁnd the eigenkets of the
many-body Hamiltonian
ˆH |Ψµ(cid:105) = Eµ |Ψµ(cid:105) ,
ˆH = ˆZ + ˆV + ...
(2.14)
14
where ˆZ, ˆV are the one-body and two-body pieces of the Hamiltonian, and in general these
terms go up to A-body interactions for a system of A particles. The many-body state is an
element of the A-body Hilbert space
|Ψ(cid:105) ∈ HA.
(2.15)
To express the many-body state in terms of single-particle quantum mechanics, let’s ﬁrst try
the A-body Hilbert space as the tensor product of A single-particle Hilbert spaces
(cid:124)
HA = H ⊗ H ⊗ ··· ⊗ H
,
(2.16)
(cid:123)(cid:122)
A
(cid:125)
where there ith single particle state is |ψi(cid:105) ∈ H. This accurately represents a many-body
state where none of the particles are interacting with each other
|Ψ(cid:105) = |ψp1ψp2 . . . ψpA(cid:105)
def
= |ψp1(cid:105) ⊗ |ψp2(cid:105) ⊗ ··· ⊗ |ψpA(cid:105) ,
(2.17)
where the positions of the kets matters, as |ψa(cid:105)|ψb(cid:105) means that particle 1 is in state |ψa(cid:105)
and particle 2 is in state |ψb(cid:105), and is in general diﬀerent from |ψb(cid:105)|ψa(cid:105). The many-body
state written as a many-body wavefunction:
(cid:104)x1, x2, . . . , xA|Ψ(cid:105) = Ψ(x1, x2, . . . , xA) = (cid:104)x1, x2, . . . , xA|ψp1ψp2 . . . ψpA(cid:105)
= ψp1(x1)ψp2(x2) . . . ψpA(xA),
(2.18)
(2.19)
where x1, x2, . . . xA represent the coordinates of the degrees of freedom (like position and
spin for example) of each particle. These are called “product states” as the wave functions
15
simply multiply together, and they form a complete A-body basis
Ψ(x1, ..., xA) =
dp1...pAψp1(x1)...ψpA(xA),
(2.20)
(cid:88)
p1,...,pA
dp1...pA = (cid:104)ψp1...ψpA|Ψ(cid:105) ,
(2.21)
where dp1...pA deﬁnes the overlap between the states.
In the case of identical particles,
interchanging any two particles in a state should leave any observable unchanged. The
particle permutation operator Pij takes two particles, i in state |ψa(cid:105) and j in state |ψb(cid:105) and
swaps them so that they occupy each other’s state. Two states are physically equivalent if
they only diﬀer by a complex phase, so
P12 |ψa(cid:105)|ψb(cid:105) = ±|ψb(cid:105)|ψa(cid:105) ,
(2.22)
gives two classes of indentical particles. Particles that are symmetric under particle inter-
change are called bosons, and particles that are antisymmetric under particle interchange
are called fermions. We will primarily be working with systems of fermions, so we intro-
duce the antisymmetrization operator because the product wave functions do not guarantee
antisymmetry
(cid:88)
ˆQ∈SA
ˆA =
1
A!
(−1)R ˆQ,
(2.23)
where A is the number of particles, SA is the symmetric group, and Q is a permutation
operator in the symmetric group, with (−1)R the associated phase of the permutation. This
will be talked about further, but for each pair of particles that are interchanged, a minus
16
sign is incurred. So if there are an even number of swaps, the sign is +, and if there are an
odd number of swaps, the sign is −. The antisymmetrizer is hermitian
and idempotent
ˆA† = ˆA,
ˆA2 = ˆA.
(2.24)
(2.25)
The antisymmetrizer projects any many-body wave function into an antisymmetric subspace.
A fermionic state is already antisymmetric, so the antisymmetrizer will act as the identity
operator when acting on a fermion wave function
|Ψ(cid:105)fermionic = ˆA|Ψ(cid:105)fermionic .
(2.26)
We can write our many-fermion ket state as the antisymmetric projection of the product ket
state space
|Ψ(cid:105)fermionic ≡ |Ψ(cid:105) =
(cid:88)
p1,...,pA
dp1...pA
ˆA|ψp1...ψpA(cid:105) .
(2.27)
Multiplying and dividing by √A!, we can rewrite this expression as
|Ψ(cid:105) =
and we will deﬁne
(2.28)
(2.29)
(cid:88)
p1,...,pA
1
√A!
dp1...pA
√A! ˆA|ψp1...ψpA(cid:105) ,
Dp1...pA ≡
1
√A!
dp1...pA,
17
and
|Φp1...pA(cid:105) ≡
√A! ˆA|ψp1...ψpA(cid:105) .
(2.30)
These many-body kets are explicitly antisymmetric, and they will form the basis for our
many-fermion basis.
If this expression is reorganized slightly, one can see that it can be written as the de-
terminant of a matrix, where the all of the permutations of diﬀerent particles in diﬀerent
single-particle states are the matrix entries, that is
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)
|Φp1...pA(cid:105) →
1
√A!
ψp1(x1)
...
. . . ψp1(xA)
. . .
...
ψpA(x1)
. . . ψpA(xA)
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) .
(2.31)
Since the determinant of a matrix is unchanged up to a sign under row/column permutation,
this representation encodes the fermionic nature of the many-body state. This is a Slater
determinant [29], and they form a complete, orthogonal and antisymmetric many-body basis
to work with. For a simple example, we can look at the two-fermion case. We will try to
reserve capital phi (Φ) as a variable representing a Slater determinant throughout this text
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)ψp1(x1) ψp1(x2)
ψp2(x1) ψp2(x2)
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)
(cid:19)
Φp1p2(x1, x2) =
1
√2
(cid:18)
=
1
√2
ψp1(x1)ψp2(x2) − ψp1(x2)ψp2(x1)
.
(2.32)
(2.33)
We see that both the antisymmetry and the Pauli exclusion principle are baked into these
18
Slater determinants
|Φp1p2...pA(cid:105) = −|Φp2p1...pA(cid:105) ,
(2.34)
and if p1 = p2 the Slater determinant is zero. So again, any fermionic wave function can be
expanded in this Slater determinant basis
(cid:88)
p1...pA
|Ψ(cid:105) =
Dp1...pA |Φp1...pA(cid:105) .
(2.35)
Looking at the coeﬃcients we see that
(cid:122)(cid:125)(cid:124)(cid:123)
antisymmetric
Hermitian
(cid:122)(cid:125)(cid:124)(cid:123)
ˆA Ψ(cid:105)
dp1...pA = (cid:104)ψp1...ψpA|
Ψ (cid:105) = (cid:104)ψp1...ψpA|
= (cid:104) ˆA(ψp1...ψpA)|Ψ(cid:105) =
1
√A! (cid:104)Φp1...pA|Ψ(cid:105) .
(2.36)
(2.37)
Thus a sign change in Φ will imply the same change in D. Also, since the Slater determinant
is zero if any two particles occupy the same state
(cid:88)
(cid:88)
|Ψ(cid:105) =
p1 = p2 =...= pA
Dp1...pA |Φp1...pA(cid:105) = A!
p1 F , while
the labels p, q, . . . represent any possible single-particle state.
The focus of the work is on inﬁnite systems, where the one-body part of the Hamiltonian
is given by the kinetic energy operator only. In second quantization it is deﬁned as
(cid:88)
pq (cid:104)p|ˆt|q(cid:105)a†paq,
ˆH0 = ˆT =
(2.114)
where the matrix elements (cid:104)p|ˆt|q(cid:105) represent the expectation value of the kinetic energy op-
erator (see the discussion below as well). The two-body interaction reads
(cid:88)
pqrs(cid:104)pq|ˆv|rs(cid:105)ASa†pa†qasar,
ˆHI = ˆV =
1
4
(2.115)
where we have deﬁned the anti-symmetrized matrix elements
(cid:104)pq|ˆv|rs(cid:105)AS = (cid:104)pq|ˆv|rs(cid:105) − (cid:104)pq|ˆv|sr(cid:105).
(2.116)
We can also deﬁne a three-body operator
ˆW =
1
36
(cid:104)pqr| ˆw|stu(cid:105)ASa†pa†qa†rauatas,
(2.117)
(cid:88)
pqrstu
47
with the anti-symmetrized matrix element
(cid:104)pqr| ˆw|stu(cid:105)AS = (cid:104)pqr| ˆw|stu(cid:105) + (cid:104)pqr| ˆw|tus(cid:105) + (cid:104)pqr| ˆw|ust(cid:105)
− (cid:104)pqr| ˆw|sut(cid:105) − (cid:104)pqr| ˆw|tsu(cid:105) − (cid:104)pqr| ˆw|uts(cid:105).
In this and the forthcoming chapters we will limit ourselves to two-body interactions at most.
Throughout this chapter and the subsequent three we will drop the subscript AS and use
only anti-symmetrized matrix elements.
Using the ansatz for the ground state |Φ0(cid:105) as new reference vacuum state, we need to
redeﬁne the anticommutation relations to
(cid:110)
a†p, aq
(cid:111)
= δpq, p, q ≤ F,
(2.118)
= δpq, p, q > F.
(2.119)
and
It is easy to see that
(cid:110)
ap, a†q
(cid:111)
ai|Φ0(cid:105) = |Φi(cid:105) (cid:54)= 0,
a†a|Φ0(cid:105) = |Φa(cid:105) (cid:54)= 0,
and
a†i|Φ0(cid:105) = 0
aa|Φ0(cid:105) = 0.
With the new reference vacuum state the Hamiltonian can be rewritten as,
(2.120)
(2.121)
ˆH = ERef + ˆHN ,
(2.122)
48
with the reference energy deﬁned as the expectation value of the Hamiltonian using the
reference state Φ0
ERef = (cid:104)Φ0| ˆH|Φ0(cid:105) =
(cid:88)
i≤F
(cid:104)i|ˆh0|i(cid:105) +
1
2
(cid:88)
ij≤F
(cid:104)ij|ˆv|ij(cid:105),
(2.123)
and the new normal-ordered Hamiltonian (all creation operators to the left of the annihilation
operators) is deﬁned as
(cid:110)
(cid:88)
pq (cid:104)p|ˆh0|q(cid:105)
a†paq
(cid:111)
(cid:88)
pqrs(cid:104)pq|ˆv|rs(cid:105)
+
1
4
(cid:110)
a†pa†qasar
(cid:111)
+
(cid:88)
pq,i≤F
(cid:104)pi|ˆv|qi(cid:105)
(cid:110)
a†paq
(cid:111)
, (2.124)
ˆHN =
where the curly brackets represent normal-ordering with respect to the new reference vacuum
state. The normal-ordered Hamiltonian can be rewritten in terms of a new one-body operator
and a two-body operator as
with
where
ˆHN = ˆFN + ˆVN ,
(cid:88)
pq (cid:104)p| ˆf|q(cid:105)
(cid:110)
a†paq
(cid:111)
,
ˆFN =
(cid:104)p| ˆf|q(cid:105) = (cid:104)p|ˆh0|q(cid:105) +
(cid:88)
i≤F
(cid:104)pi|ˆv|qi(cid:105).
(2.125)
(2.126)
(2.127)
The last term on the right hand side represents a medium modiﬁcation to the single-particle
Hamiltonian due to the two-body interaction. Finally, the two-body interaction is given by
(cid:110)
(cid:88)
pqrs(cid:104)pq|ˆv|rs(cid:105)
ˆVN =
1
4
(cid:111)
a†pa†qasar
.
49
2.12 Conﬁguration Interaction
With this in place, let’s look again at the many-body Schr¨odinger equation.
ˆH |Ψ(cid:105) = E |Ψ(cid:105)
(2.128)
and we expand the solution to the Schr¨odinger equation in the basis of our complete and
orthonormal antisymmetrized Slater determinants
(cid:88)
(cid:88)
i
i∈all SD’s
(cid:88)
i
ˆH
|Ψ(cid:105) =
Ci |Φi(cid:105) ,
Ci |Φi(cid:105) = E
Ci |Φi(cid:105) .
(2.129)
(2.130)
As stated in the Slater rules section, we can project Eqn. (2.128) onto (cid:104)Φj|, yielding
(cid:88)
(cid:88)
i
(cid:88)
i
(cid:104)Φj| ˆH
Ci |Φi(cid:105) = (cid:104)Φj| E
HijCi = ECi,
Ci |Φi(cid:105) ,
(2.131)
(2.132)
i
which is exactly the index notation for ←→H (cid:126)C = E (cid:126)C, which is the classic statement of the
eigenvalue problem in linear algebra, where the Hamiltonian can be written as a matrix with
matrix elements
Hij = (cid:104)Φi| ˆH|Φj(cid:105) .
To begin solving this equation, ﬁrst a single-particle basis truncation must be made. Since
the many-body basis is built up from a single-particle basis, this ensures there are a ﬁnite
50
number of Slater determinants, and thus a ﬁnite number of matrix elements. Once this
matrix is ﬁnite, the eigenvalues and eigenvectors can be found by diagonalizing this matrix,
giving the exact solution within this truncated space. Recomputing this solution for larger
and larger single-particle basis set cutoﬀs will create a series of solutions that can hopefully
generate a smooth curve and an inﬁnite basis limit can be extrapolated to. Thus the many-
body problem is solved.
Well, not quite. Since the size of the Hamiltonian matrix is N × N where N is the
number of Slater determinants, and since the Slater determinants are generated by creating
all permutations in the symmetry group Sn for n single-particle states which grows factorially,
this matrix gets very large, very fast. So computing these matrix diagonalizations with even
the largest modern supercomputers quickly becomes impossible even for systems on the
order of 10 particles. This means that further approximations are necessary to compute
properties of larger quantum systems. Along with the truncation to the single-particle basis,
we can truncate the basis of Slater determinants. While it is valid to take the list of Slater
determinants and throw a large fraction away so that the Hamiltonian matrix is smaller, not
all truncations are equal. This is where the machinery developed in the Fermi level formalism
section comes into use. In most many-body problems, not all conﬁgurations of particles and
states are equally important. For example, the probability of ﬁnding all A particles in the
A highest lying states in your basis is usually vanishingly small. So the scheme developed
in terms of particle and hole conﬁgurations with respect to some reference state keep the
“important” states at the front of our attention. Then, if truncations need to be made, we
can hopefully truncate Slater determinants which do not have much overlap with the ground
state. This of course assumes that an “important” reference state can indeed be picked out,
and this procedure will be expanded on further in the Hartree-Fock section.
51
(cid:88)
i,a
(cid:88)
i,j,a,b
(cid:88)
C
i1,i2,...,iA
a1,a2,...aA
So the complete expansion of a many body state as in Eqn. (2.129) can be rewritten as
|ΨF CI(cid:105) = C0 |Φ0(cid:105) +
i |Φa
Ca
i (cid:105) +
ij |Φab
Cab
ij (cid:105) + ··· +
a1a2...aA
i1i2...iA |Φ
a1a2...aA
i1i2...iA (cid:105) ,
(2.133)
where F CI stands for full conﬁguration interaction, meaning that the full spaces of inter-
acting conﬁgurations is being used in this expansion. This is a complete A-body basis, just
reorganizing the Slater determinants in terms of how many excitations away from the refer-
ence state they are. The ﬁrst approximation is a ﬁnite cutoﬀ for the inﬁnite single-particle
basis. Once the single-particle basis has been truncated, the above expression gives the exact
answer in this subspace, as FCI includes all Slater determinants up to A-body excitations,
at which point the series naturally truncates. The next natural approximation, is to not
include every Slater determinant. For example, if only singles and doubles excitations are
included, then we get the conﬁguration interaction singles doubles (CISD) approximation,
which looks like:
|ΨCISD(cid:105) = C0 |Φ0(cid:105) +
(cid:88)
i,a
Ca
i |Φa
i (cid:105) +
(cid:88)
i,j,a,b
Cab
ij |Φab
ij (cid:105) .
(2.134)
While this is a conceptually nice place to truncate the series, it turns to be a pretty poor
way to include correlations into the target many-body state for a given computational cost.
The various ways to optimize these truncated many-body correlations is a widely studied
topic in its own right. Unfortunately, there is no one many-body method to rule them all,
as each approximation has its own strengths and weaknesses. This spans from FCI which
includes every many-body correlation, down to Hartree-Fock which includes only a single
52
Slater determinant.
2.13 Hartree-Fock Theory
Hartree-Fock (HF) theory [35, 36], is an algorithm for ﬁnding an approximative expression
for the ground state of a given Hamiltonian. The basic ingredients contain a single-particle
basis {ψα} deﬁned by the solution of the following eigenvalue problem
ˆhHFψα = εαψα,
(2.135)
with the Hartree-Fock Hamiltonian deﬁned as
ˆhHF = ˆt + ˆuext + ˆuHF.
(2.136)
The term ˆuHF is a single-particle potential to be determined by the HF algorithm. The
HF algorithm means to select ˆuHF in order to have
(cid:104) ˆH(cid:105) = EHF = (cid:104)ΦHF
0
| ˆH|ΦHF
0
(cid:105),
(2.137)
as a local minimum with a Slater determinant ΦHF
0
being the ansatz for the ground state.
The variational principle ensures that EHF ≥ E0, with E0 representing the exact ground
state energy.
We will show that the Hartree-Fock Hamiltonian ˆhHF equals our deﬁnition of the operator
ˆf discussed in connection with the new deﬁnition of the normal-ordered Hamiltonian, that
53
is we have, for a speciﬁc matrix element
(cid:104)p|ˆhHF|q(cid:105) = (cid:104)p| ˆf|q(cid:105) = (cid:104)p|ˆt + ˆuext|q(cid:105) +
(cid:88)
i≤F
(cid:104)pi| ˆV |qi(cid:105),
(2.138)
meaning that
(cid:88)
i≤F
(cid:104)pi| ˆV |qi(cid:105).
(cid:104)p|ˆuHF|q(cid:105) =
(2.139)
The so-called Hartree-Fock potential ˆuHF adds an explicit medium dependence due to the
summation over all single-particle states below the Fermi level F . It brings also in an explicit
dependence on the two-body interaction (in nuclear physics we can also have complicated
three- or higher-body forces). The two-body interaction, with its contribution from the
other bystanding fermions, creates an eﬀective mean ﬁeld in which a given fermion moves,
in addition to the external potential ˆuext which conﬁnes the motion of the fermion. For
systems like nuclei or inﬁnite nuclear matter, there is no external conﬁning potential. Nuclei
and nuclear matter are examples of self-bound systems, where the binding arises due to the
intrinsic nature of the strong force. For nuclear systems thus, there would be no external
one-body potential in the Hartree-Fock Hamiltonian.
Another possibility is to expand the single-particle functions in a known basis and vary
the coeﬃcients, that is, the new single-particle wave function is written as a linear expansion
in terms of a ﬁxed chosen orthogonal basis (for example the well-known harmonic oscillator
functions or the hydrogen-like functions, etc). We deﬁne our new Hartree-Fock single-particle
basis by performing a unitary transformation on our previous basis (labelled with Greek
indices) as
ψHF
p =
(cid:88)
λ
54
Cpλφλ.
(2.140)
In this case we vary the coeﬃcients Cpλ. If the basis has inﬁnitely many terms, we need to
truncate the above sum. We assume that the basis φλ is orthogonal. A unitary transforma-
tion keeps the orthogonality, which is desired.
It is normal to choose a single-particle basis deﬁned as the eigenfunctions of parts of the
full Hamiltonian. The typical situation consists of the solutions of the one-body part of the
Hamiltonian, that is we have
ˆh0φλ = λφλ.
(2.141)
For inﬁnite nuclear matter, ˆh0 is given by the kinetic energy operator and the states are given
by plane wave functions. Due to the translational invariance of the two-body interaction,
the Hartree-Fock single-particle eigenstates are also given by the same functions. Thus, for
inﬁnite matter it is only the single-particle energies that change when we solve the Hartree-
Fock equations.
The single-particle wave functions φλ(r), deﬁned by the quantum numbers λ and r are
deﬁned as the overlap
φλ(r) = (cid:104)r|λ(cid:105).
(2.142)
In our discussions we will use our deﬁnitions of single-particle states above and below the
Fermi (F ).
We use Greek letters to refer to our original single-particle basis. The expectation value
for the energy with the ansatz Φ0 for the ground state reads
(cid:88)
µ≤F
(cid:88)
µ,ν≤F
E[Φ0] =
(cid:104)µ|h|µ(cid:105) +
1
2
(cid:104)µν|ˆv|µν(cid:105).
(2.143)
Now we are interested in deﬁning a new basis deﬁned in terms of a chosen basis as deﬁned
55
in Eq. (2.140). We deﬁne the energy functional as
E[ΦHF ] =
(cid:104)i|h|i(cid:105) +
1
2
(cid:88)
i≤F
(cid:88)
ij≤F
(cid:104)ij|ˆv|ij(cid:105),
(2.144)
where ΦHF is the new Slater determinant deﬁned by the new basis of Eq. (2.140).
Using Eq. (2.140) we can rewrite Eq. (2.144) as
(cid:88)
(cid:88)
i≤F
αβ
(cid:88)
(cid:88)
ij≤F
αβγδ
E[Ψ] =
C∗iαCiβ(cid:104)α|h|β(cid:105) +
1
2
C∗iαC∗jβCiγCjδ(cid:104)αβ|ˆv|γδ(cid:105).
(2.145)
In order to ﬁnd the variational minimum of the above functional, we introduce a set of
Lagrange multipliers, noting that since (cid:104)i|j(cid:105) = δi,j and (cid:104)α|β(cid:105) = δα,β, the coeﬃcients Ciγ
obey the relation
(cid:104)i|j(cid:105) = δi,j =
C∗iαCiβ(cid:104)α|β(cid:105) =
C∗iαCiα,
(2.146)
(cid:88)
αβ
(cid:88)
α
(cid:88)
(cid:88)
i≤F
α
which allows us to deﬁne a functional to be minimized that reads
F [ΦHF ] = E[ΦHF ] −
i
C∗iαCiα.
(2.147)
Minimizing with respect to C∗iα (the equations for C∗iα and Ciα can be written as two
independent equations) we obtain
E[ΦHF ] −
= 0,
C∗jαCjα
(cid:88)
j
j
(cid:88)
α
(2.148)
d
dC∗iα
which yields for every single-particle state i and index α (recalling that the coeﬃcients Ciα
are matrix elements of a unitary matrix, or orthogonal for a real symmetric matrix) the
56
following Hartree-Fock equations
(cid:88)
β
Ciβ(cid:104)α|h|β(cid:105) +
(cid:88)
(cid:88)
j≤F
βγδ
C∗jβCjδCiγ(cid:104)αβ|ˆv|γδ(cid:105) = HF
i Ciα.
(2.149)
We can rewrite this equation as (changing dummy variables)
(cid:104)α|h|β(cid:105) +
(cid:88)
β
(cid:88)
(cid:88)
j≤F
γδ
Ciβ = HF
C∗jγCjδ(cid:104)αγ|ˆv|βδ(cid:105)
i Ciα.
(2.150)
Note that the sums over Greek indices run over the number of basis set functions (in principle
an inﬁnite number).
Deﬁning
hHF
αβ = (cid:104)α|h|β(cid:105) +
(cid:88)
(cid:88)
j≤F
γδ
C∗jγCjδ(cid:104)αγ|ˆv|βδ(cid:105),
(2.151)
we can rewrite the new equations as
(cid:88)
β
hHF
αβ Ciβ = HF
i Ciα.
(2.152)
The latter is nothing but a standard eigenvalue problem. Our Hartree-Fock matrix is thus
ˆhHF
αβ = (cid:104)α|ˆh0|β(cid:105) +
C∗jγCjδ(cid:104)αγ|ˆv|βδ(cid:105).
(2.153)
(cid:88)
(cid:88)
j≤F
γδ
The Hartree-Fock equations are solved in an iterative way starting with a guess for the
coeﬃcients Cjγ = δj,γ and solving the equations by diagonalization till the new single-particle
energies HF
i
do not change anymore by a user-deﬁned small quantity.
Normally we assume that the single-particle basis |β(cid:105) forms an eigenbasis for the operator
57
ˆh0, meaning that the Hartree-Fock matrix becomes
(cid:88)
(cid:88)
j≤F
γδ
ˆhHF
αβ = αδα,β +
C∗jγCjδ(cid:104)αγ|ˆv|βδ(cid:105).
(2.154)
2.14 Many-Body Perturbation Theory
Hartree-Fock theory is incredibly powerful considering how simple the idea is: transform
the single-particle basis to optimize a single Slater determinant. This optimization problem
can be thought of as summing the interactions from the surrounding particles to create a
sort of external “mean ﬁeld” that each particle feels. In some cases, this mean ﬁeld approach
is suﬃcient to answer the physics questions being asked. However, the heart of many-body
theory is the correlation between particles that are purely many-body in nature, meaning
that they cannot be described in the independent particle picture. Historically, many-body
perturbation theory (MBPT) [37, 38, 33, 16] has been a ﬁrst attempt to build these many-
body correlations on top of a reference state. Like perturbation theory in other branches of
physics, or in single particle quantum mechanics, it is assumed that the problem at hand
is a small “perturbation” away from a reference problem. In the many-body case, starting
from a good reference state is essential, as a bad starting point will require many additional
corrections, or in some cases adding additional corrections will not work at all. To set the
stage for deriving MBPT, we assume here that we are only interested in the non-degenerate
ground state of a given system and expand the exact wave function in terms of a series of
Slater determinants
∞(cid:88)
m=1
Cm|Φm(cid:105),
|Ψ0(cid:105) = |Φ0(cid:105) +
58
(2.155)
where we have assumed that the true ground state is dominated by the solution of the
unperturbed problem, that is
ˆH0|Φ0(cid:105) = W0|Φ0(cid:105),
and that the full Hamiltonian is given by this term plus a small interaction term
ˆH = ˆH0 + ˆHI .
(2.156)
(2.157)
The state |Ψ0(cid:105) is not normalized and we employ again intermediate normalization via
(cid:104)Φ0|Ψ0(cid:105) = 1.
The Schr¨odinger equation is given by
ˆH|Ψ0(cid:105) = E|Ψ0(cid:105),
and multiplying the latter from the left with (cid:104)Φ0| gives
(cid:104)Φ0| ˆH|Ψ0(cid:105) = E(cid:104)Φ0|Ψ0(cid:105) = E,
and subtracting from this equation
(2.158)
(2.159)
(cid:104)Ψ0| ˆH0|Φ0(cid:105) = W0(cid:104)Ψ0|Φ0(cid:105) = W0,
(2.160)
and using the fact that the operators ˆH and ˆH0 are hermitian results in
∆E = E − W0 = (cid:104)Φ0| ˆHI|Ψ0(cid:105),
(2.161)
59
which is an exact result. The total energy can be separated into two terms
E = ERef + ∆E,
(2.162)
where ∆E is the correlation energy, and the reference energy is given by
ERef = (cid:104)Φ0| ˆH|Φ0(cid:105).
(2.163)
Equation (2.161) forms the starting point for all perturbative derivations. However, as
it stands it represents nothing but a mere formal rewriting of Schr¨odinger’s equation and is
not of much practical use. The exact wave function |Ψ0(cid:105) is unknown. In order to obtain a
perturbative expansion, we need to expand the exact wave function in terms of the interaction
ˆHI .
Here we have assumed that our model space deﬁned by the operator ˆP is one-dimensional,
meaning that
and
ˆP = |Φ0(cid:105)(cid:104)Φ0|,
∞(cid:88)
ˆQ =
|Φm(cid:105)(cid:104)Φm|.
m=1
We can thus rewrite the exact wave function as
|Ψ0(cid:105) = ( ˆP + ˆQ)|Ψ0(cid:105) = |Φ0(cid:105) + ˆQ|Ψ0(cid:105).
(2.164)
(2.165)
(2.166)
Going back to the Schr¨odinger equation, we can rewrite it, adding and a subtracting a term
60
ω|Ψ0(cid:105) as
(cid:17)
(cid:16)
ω − ˆH0
|Ψ0(cid:105) =
(cid:16)
ω − E + ˆHI
(cid:17)
|Ψ0(cid:105),
(2.167)
where ω is an energy variable to be speciﬁed later.
We assume also that the resolvent of
ω − ˆH0
deﬁnes the unperturbed Green’s function as
(cid:17)
exits, that is it has an inverse which
(cid:16)
(cid:17)−1
(cid:16)
ω − ˆH0
|Ψ0(cid:105) =
1
ω − ˆH0
(cid:16)
=
(cid:17).
1
ω − ˆH0
(cid:17)
|Ψ0(cid:105),
ω − E + ˆHI
(cid:16)
We can rewrite Schr¨odinger’s equation as
(2.168)
(2.169)
(2.170)
and multiplying from the left with ˆQ results in
(cid:16)
ω − E + ˆHI
(cid:17)
|Ψ0(cid:105),
ˆQ|Ψ0(cid:105) =
ˆQ
ω − ˆH0
which is possible since we have deﬁned the operator ˆQ in terms of the eigenfunctions of ˆH0.
Since these operators commute we have
(cid:16)
ˆQ
(cid:17) ˆQ = ˆQ
(cid:16)
(cid:17) =
(cid:17).
ˆQ(cid:16)
ω − ˆH0
1
ω − ˆH0
1
ω − ˆH0
(2.171)
With these deﬁnitions we can in turn deﬁne the wave function as
(cid:16)
(cid:17)
ω − E + ˆHI
|Ψ0(cid:105).
(2.172)
|Ψ0(cid:105) = |Φ0(cid:105) +
ˆQ
ω − ˆH0
61
So far, this is just a reorganization of the Schr¨odinger equation. It is a non-linear equation
in two unknown quantities, the energy E and the exact wave function |Ψ0(cid:105). We can however
start with a guess for |Ψ0(cid:105) on the right hand side of the last equation.
The most common choice is to start with the function which is expected to exhibit the
largest overlap with the wave function we are searching after, namely |Φ0(cid:105). This can again
be inserted in the solution for |Ψ0(cid:105) in an iterative fashion and if we continue along these
lines we end up with
(cid:40) ˆQ
∞(cid:88)
(cid:16)
(cid:17)(cid:41)i
|Ψ0(cid:105) =
i=0
ω − ˆH0
ω − E + ˆHI
|Φ0(cid:105),
(2.173)
for the wave function and
∆E =
∞(cid:88)
i=0
(cid:104)Φ0| ˆHI
(cid:40) ˆQ
(cid:16)
ω − ˆH0
(cid:17)(cid:41)i
ω − E + ˆHI
|Φ0(cid:105),
(2.174)
which is now a perturbative expansion of the exact energy in terms of the interaction ˆHI
and the unperturbed wave function |Ψ0(cid:105).
In our equations for |Ψ0(cid:105) and ∆E in terms of the unperturbed solutions |Φi(cid:105) we have
still an undetermined parameter ω and a dependency on the exact energy E. Not much has
been gained thus from a practical computational point of view.
In Brillouin-Wigner perturbation theory [16, 39] it is customary to set ω = E. This
62
results in the following perturbative expansion for the energy ∆E
∞(cid:88)
(cid:32)
(cid:104)Φ0| ˆHI
i=0
(cid:40) ˆQ
(cid:16)
ω − ˆH0
ˆQ
∆E =
(cid:17)(cid:41)i
ω − E + ˆHI
|Φ0(cid:105) =
(cid:104)Φ0|
ˆHI + ˆHI
ˆHI + ˆHI
ˆQ
E − ˆH0
ˆHI
ˆQ
E − ˆH0
E − ˆH0
(cid:33)
ˆHI + . . .
|Φ0(cid:105).
(2.175)
(2.176)
This expression depends however on the exact energy E and is again not very convenient
from a practical point of view. It can obviously be solved iteratively, by starting with a guess
for E and then solve till some kind of self-consistency criterion has been reached.
Deﬁning e = E − ˆH0 and recalling that ˆH0 commutes with ˆQ by construction and that ˆQ
is an idempotent operator ˆQ2 = ˆQ, we can rewrite the denominator in the above expansion
for ∆E as
ˆQ
1
ˆe − ˆQ ˆHI ˆQ
= ˆQ
(cid:20)1
ˆe
+
1
ˆe
ˆQ ˆHI ˆQ
1
ˆe
+
1
ˆe
ˆQ ˆHI ˆQ
1
ˆe
ˆQ ˆHI ˆQ
1
ˆe
+ . . .
(cid:21)
ˆQ.
Inserted in the expression for ∆E, we obtain
∆E = (cid:104)Φ0| ˆHI + ˆHI ˆQ
1
E − ˆH0 − ˆQ ˆHI ˆQ
ˆQ ˆHI|Φ0(cid:105).
(2.177)
(2.178)
In Rayleigh-Schr¨odinger (RS) perturbation theory [40, 41, 16] we set ω = W0 and obtain the
following expression for the energy diﬀerence
∞(cid:88)
i=0
∆E =
(cid:40)
(cid:32)
(cid:104)Φ0| ˆHI
ˆQ
W0 − ˆH0
ˆQ
(cid:17)(cid:41)i
(cid:16)
ˆHI − ∆E
|Φ0(cid:105)
(cid:33)
=(cid:104)Φ0|
ˆHI + ˆHI
( ˆHI − ∆E) + . . .
|Φ0(cid:105).
W0 − ˆH0
63
(2.179)
(2.180)
The operator ˆQ commutes with ˆH0 and since ∆E is a constant we obtain that
ˆQ∆E|Φ0(cid:105) = ˆQ∆E| ˆQΦ0(cid:105) = 0.
Inserting this result in the expression for the energy gives us
(cid:32)
∆E = (cid:104)Φ0|
ˆHI + ˆHI
ˆQ
W0 − ˆH0
ˆHI + ˆHI
ˆQ
W0 − ˆH0
( ˆHI − ∆E)
ˆQ
W0 − ˆH0
ˆHI + . . .
(2.181)
(cid:33)
|Φ0(cid:105).
(2.182)
We can now perturbatively expand this expression in terms of the interaction ˆHI , which
is assumed to be small. We obtain then
∆E =
with the following expression for ∆E(i)
∞(cid:88)
i=1
∆E(i),
∆E(1) = (cid:104)Φ0| ˆHI|Φ0(cid:105),
which is just the contribution to ﬁrst order in perturbation theory,
∆E(2) = (cid:104)Φ0| ˆHI
ˆQ
W0 − ˆH0
ˆHI|Φ0(cid:105),
(2.183)
(2.184)
(2.185)
which is the contribution to second order. There exists a formal theory for the calculation of
each additional term, see for example Ref. [16], where a diagrammatic method is described to
generate any order of MBPT. Inserting in the ˆQ space operator and the energy denominators
in Eqn. (2.185) we get the expression for MBPT(2), the energy correction for many-body
64
perturbation theory to second order,
∆E(2) =
1
22
(cid:88)
(cid:88)
ij≤F
ab>F
(cid:104)ij|ˆv|ab(cid:105)(cid:104)ab|ˆv|ij(cid:105)
εi + εj − εa − εb
.
(2.186)
In the expressions for the various diagrams the quantity ε denotes the single-particle
energies deﬁned by H0. Many-body perturbation theory is quite powerful, and can provide
corrections quite accurately for many systems. Unfortunately, most nuclear physics appli-
cations do not fall in this regime where MBPT is accurate due to that large short range
correlations of the nucleon-nucelon potential. Any system for which the interactions are too
large compared to the non-interacting mean-ﬁeld model will not be meaningfully captured
by MBPT even at high orders of correction. This will be illustrated in Chapter 3, where the
simple pairing model is examined.
2.15
In-Medium Similarity Renormalization Group
So far we have covered three important methods in quantum many-body theory. First,
was Full conﬁguration interaction (FCI), where the only approximation is the necessary
truncation to the single-particle basis. The many-body basis is then constructed from the
single-particle basis, but FCI does not make any additional truncations to the many-body
space. Next, was Hartree-Fock (HF) mean ﬁeld theory, which optimizes the many-body
ground state energy in a single Slater determinant by performing a unitary transformation
on the single-particle basis. Last, was many-body perturbation theory (MBPT), which exists
between a simple mean ﬁeld model and the full many-body solution. This is the regime
where the vast majority of many-body physics is done, as many interesting problems are
65
too computationally expensive for FCI, but too correlated for Hartree-Fock to be suﬃcient.
As a result, many diﬀerent approximations to the many-body problem have been developed
with various strengths and weaknesses to target diﬀerent applications. The strength of ab
initio many-body methods is that by maintaining the essential degrees of freedom, all of the
approximations can be extended to recover the full solution. Of course the trade-oﬀ is then
that these extensions towards the exact solutions again become prohibitively expensive, but
this allows for a choice in the trade-oﬀ between accuracy in computational cost. Additionally,
in nuclear physics the nuclear potentials can be developed in a similar way that allows for
increased accuracy (at increasing computational cost), but this work will focus on various
approximations to the many-body methods rather than the input potentials. In particular,
most of the calculations in this work are done with an ab initio many-body method called
coupled cluster (CC) theory, which is explained in detail in Chapter 4. Another ab initio
many-body method similar to coupled cluster is in-medium similarity renormalization group
(IM-SRG) [42, 43, 44, 45, 46, 47, 48, 3].
The renormalization group is a tool that has been used in physics for many decades,
which allows physical quantities of interest to be examined at diﬀerent distance or energy
scales and has been essential in the development of quantum electrodynamics and quantum
chromodynamics. In nuclear physics most realistic nuclear potentials have a sharp repulsive
core which can lead to divergences in calculating matrix elements, generating strong oﬀ-
diagonal contributions as low momentum modes are coupled to high momentum modes.
However for certain physical quantities, like the ground state energy of an atomic nucleus,
the low energy of the system indicates that the nucleons should not probe the very short
range distance scales of this repulsive core. Similarity renormalization group (SRG) has
had success in “softening” the repulsive of the nuclear potential, by driving the momenum
66
space interaction matrix to a band diagonal form, decoupling the high momentum from low
momentum modes while maintaining accuracy of the target observables [4, 5].
In-medium similarity renormalization group (IM-SRG) takes this philosophy of decou-
pling distance scales and applies it to the “medium” of a particular reference state. The idea
is that for a matrix problem in the many-body basis, transforming the matrix to a form where
the ground state is decoupled from the rest of the matrix will give the ground state eigen-
value. In the case of IM-SRG, the Hamiltonian is normal ordered with respect to a reference
state and the ground state energy is isolated by a continuous unitary transformation.
A unitary transformation U is an isomorphism between two Hilbert spaces H1, H2 that
preserves the inner product,
U : H1 → H2,
(cid:104)U x, U y(cid:105)H1
= (cid:104)x, y(cid:105)H2
,
∀x, y ∈ H1. Similarly, for an antiunitary transformation,
(cid:104)U x, U y(cid:105) = (cid:104)x, y(cid:105)∗ = (cid:104)y, x(cid:105) .
Thus for a unitary transformation, in bra-ket notation,
(cid:104)U x|U y(cid:105) = (cid:104)x|U†U|y(cid:105) = (cid:104)x|y(cid:105) ,
=⇒ U†U = U U† = 1.
67
(2.187)
(2.188)
(2.189)
(2.190)
(2.191)
A continuous unitary transformation is a unitary transformation parametrized by some con-
tinuous parameter s, such that U (s)U (s)† = 1. This generates a unitarily transformed
Hamiltonian for all points s,
H(s) = U†(s)HU (s).
(2.192)
The transformation is implemented by solving a coupled set of ﬂow equations for the matrix
elements for the Hamiltonian, which we can ﬁnd by taking the derivative of Eqn. (2.192),
dH(s)
ds
=
dU†(s)
ds
HU (s) + U†(s)H
dU (s)
ds
,
and the derivative of Eqn. (2.191)
dU†(s)
ds
U (s) + U†(s)
dU (s)
ds
= 0.
From here, we write down the generator of the transformation η as
η(s) =
dU†(s)
ds
U (s) = −U†(s)
dU (s)
ds
,
which leads to the ﬂow equations as
=(cid:2)η(s), H(s)(cid:3).
dH(s)
ds
(2.193)
(2.194)
(2.195)
(2.196)
For actual calculations, an explicit expression from the transformation U (s) is rarely written
out. Instead, the generator η deﬁnes the unitary transformation. To actually implement
68
this, we partition the Hamiltonian as
H(s) = Hd(s) + Hod(s),
(2.197)
where these are the diagonal and oﬀ-diagonal components of the matrix. The evolution
with the continuous ﬂow parameter s is again The choice of the generator ﬁrst suggested by
Wegner [49],
guarantees
η(s) =(cid:2)Hd(s), H(s)(cid:3) =(cid:2)Hd(s), Hod(s)(cid:3),
(cid:18)(cid:0)Hod(cid:1)2
= 2Tr(cid:0)η2(cid:1) = −2Tr(cid:0)η†η(cid:1)
(cid:19)
≤ 0,
d
ds
Tr
(2.198)
(2.199)
which demonstrates that Hod decays with increasing s which is precisely what is needed
to decouple the high and low momentum modes[50]. Analyzing the ﬂow equations in the
eigenbasis of Hd(s) and deﬁning Hd
ii(s) ≡ i one can show that
ij (s) ∼ e−s(i−j )2
Hod
Hod
ij (0).
(2.200)
However, this can lead to stiﬀ ODE’s, so a more common generator is the White generator
[51]
which gives uniform surpression
ηij(s) =
Hod
ij (s)
i − j
,
Hod
ij (s) ∼ e−sHod
ij (0).
69
(2.201)
(2.202)
While SRG is typically used to soften nuclear potentials with a repulsive core, an alter-
native is to perform the SRG evolution in-medium (IM-SRG) for each A-body system of
interest. Starting from a general second-quantized Hamiltonian with two- and three-body
interactions
(cid:88)
qr
H =
Tqra†qar +
1
2!2
(cid:88)
qrst
(2)
qrsta†qa†ratas +
V
1
3!2
(cid:88)
qrstuv
(3)
qrstuva†qa†ra†savauat + . . .
V
(2.203)
All operators can be normal-ordered with respect to a ﬁnite-density Fermi vacuum |Φ(cid:105) (e.g.
the Hartree-Fock ground state), as opposed to the zero particle vacuum. Wick’s theorem
can then be used to exactly write H as
(cid:88)
qr
(cid:88)
qrst
H = E +
fqr{a†qar}+
1
4
Γqrst{a†qa†ratas}+
1
36
(cid:88)
qrstuv
Wqrstuv{a†qa†ra†savauat}, (2.204)
where strings of normal ordered operators obey
(cid:104)Φ|{a†q . . . ar}|Φ(cid:105) = 0,
(2.205)
and the terms in (2.204) are given by
(cid:88)
E =
Tqqnq +
q
fqr = Tqr +
(cid:88)
(cid:88)
s
u
Γqrst = V
(2)
qrst +
Wqrtsuv = V
(3)
qrstuv,
(cid:88)
qr
1
2
V
(2)
qsrsns +
V
(2)
qrqrnqnr +
(cid:88)
st
1
2
(cid:88)
qrs
1
6
(3)
qrsqrsnqnrns,
V
(2.206)
(3)
V
qstrstnsnt,
(2.207)
(2.208)
(2.209)
(3)
V
qrustunu,
70
where nq = θ(F − q) are the occupation numbers in the reference state |Φ(cid:105). Notice
that the normal ordered 0-,1-, and 2-body terms include contributions from the three-body
interaction V (3) through sums over the occupied single-particle states in the reference state
|Φ(cid:105). Neglecting the residual three-body interaction leads to the normal-ordered two-body
approximation (NO2B) which has shown to be an excellent approximation in nuclear systems.
Truncating the in-medium SRG equations to normal-ordered two-body operators is denoted
IM-SRG(2). Using this normal ordered Hamiltonian and using Wick’s theorem on Eqn.
(2.196) with H(s) = E0(s) + f (s) + Γ(s) and truncating η(s) = η(1)(s) + η(2)(s) to two-body
yields the coupled IM-SRG(2) equations
dE0
ds
=
ηqrfrq(nq − nr) +
(cid:88)
qrst
1
2
ηqrstΓstqrnqnr ¯ns¯nt,
(2.210)
dfqr
ds
=
(1 + Pqr)ηqsfsr +
(cid:88)
st
(ns − nt)(ηstΓtqsr − fstηtqsr)
(nsnt¯nu + ¯ns¯ntnu)(1 + Pqr)ηuqstΓstur,
dΓqrst
ds
=
(1 − Pqr)(ηquΓurst − fquηurst)
(cid:88)
qr
(cid:88)
(cid:88)
s
+
stu
u
(cid:88)
(cid:88)
(cid:88)
(cid:88)
u
1
2
−
+
uv
−
uv
(2.211)
(2.212)
(1 − Pst)(ηusΓqrut − fusηqrut)
(1 − nu − nv)(ηqruvΓuvst − Γqruvηuvst)
(nu − nv)(1 − Pqr)(1 − Pst)ηvrutΓuqvs,
71
where the ¯nr ≡ (1 − nr), Pqr is the permutation operator and s dependence has been made
implicit to clear up visual clutter. The White generator is then
(cid:88)
ai
η =
(cid:88)
abij
fai
fa − fi{a†aai} +
1
4
Γabij
fa + fb − fi − fj {a†aa†bajai} − H.c.,
(2.213)
where fa = faa are the Møller-Plesset energy denominators. Unfortunately, these equations
can be very sensitive, as small amounts of numerical error can break the unitarity of the
transformation. This means that to solve these equations, a high-order diﬀerential equation
solver is typically needed. These solvers need to store many copies of the solution vector to
maintain accuracy, and these copies rapidly increase the storage requirements. Fortunately,
the Magnus expansion can help out here, ensuring that unitarity is preserved at every step
in the diﬀerential equation.
2.16 The Magnus Formulation of IM-SRG
The starting point of the Magnus formulation [50] of IM-SRG is once again taking the
derivative of the unitarity condition U (s)U†(s) = U†(s)U (s) = 1,
dU (s)
ds
U†(s) = −U (s)
dU†(s)
ds
.
(2.214)
and multiply Eqn. (2.214) on the right by U (s) to yield the
Now deﬁne η ≡ U (s)
diﬀerential equation
dU†(s)
ds
dU (s)
ds
dU†(s)
U†(s)U (s) = −U (s)
=⇒
ds
= −η(s)U (s),
dU (s)
ds
U (s),
(2.215)
(2.216)
72
with the boundary condition U (0) = 1. To get some intuition for this diﬀerential equation, we
look to a familiar unitary transformation, like the time evolution operator, the Hamiltonian.
U (t) = e−iHt.
Taking the time derivative yields
dU (t)
dt
= −iHe−iHt = −iHU (t).
(2.217)
(2.218)
This is true when H is independent of t and explains why the solution is so compact. If we
look at U (s) = e−ηs, the derivative would be
dU (s)
ds
= −ηe−ηs = −ηU (s),
(2.219)
which would be a nice solution to the diﬀerential equation. So it becomes clear that the
s dependence in η(s) makes things more complicated. If we had U (s) = e−η(s), then the
derivative would be
dU (s)
ds
dη(s)
ds
= −
e−η(s).
(2.220)
Thus to get the solution we want, we need the anti-derivative of η(s) to be exponentiated;
so something like Exp(−
matrix exponential, which is really just short hand for the polynomial series. So terms like
0 η(s(cid:48))ds(cid:48)). But here, there are all sorts of issues since this is a
(cid:82) s
(cid:90) s
(cid:90) s
1
n!
0 ···
0
η(s(cid:48)1) . . . η(s(cid:48)n)ds1 . . . dsn,
(2.221)
73
arise. And here, unless all of the η terms at any value of s commute, the order matters. This
can be formally integrated as the time-ordered exponential
(cid:8)e−
0 η(s(cid:48))ds(cid:48)(cid:9)
(cid:82) s
≡ 1 −
U (s) = Ts
(cid:90) s
0
(cid:90) s
(cid:90) ds(cid:48)
0
0
ds(cid:48)η(s(cid:48)) +
ds(cid:48)
ds(cid:48)(cid:48)η(s(cid:48))η(s(cid:48)(cid:48)) + . . .
(2.222)
This is not useful from a practical point of view. The Magnus expansion [52] is the statement
that given a few technical requirements on η(s), a solution of the form
U (s) = eΩ(s)
(2.223)
exists, where Ω†(s) = −Ω(s) and Ω(0) = 0. This lines up with the previously stated
boundary condition of U (0) = 1, which is satisﬁed by Ω(0) = 0. The anti-Hermitian property
of Ω is necessary since for any unitary operator U to be expressed as exponentiated operator
Ω requires that the exponentiated operator Ω be anti-Hermitian,
U U† = 1
= eΩeΩ† = eΩ+Ω†e[Ω,Ω†].
?
(2.224)
This expression will be satisﬁed as long as Ω† = −Ω, since
(cid:2)Ω, Ω†(cid:3) = ΩΩ† − Ω†Ω = −Ω2 + Ω2 = 0
(cid:2)Ω,Ω†(cid:3)
=⇒ eΩ+Ω†e
= e0e0 = 1.
(2.225)
This is why the time evolution operator eiHt has the characteristic phase i. Since H is
Hermitian, the i is needed in the exponential to make the exponent anti-Hermitian overall
to ensure the unitarity of the transformation.
74
In previous applications of the Magnus expansion, Ω(s) is expanded in powers of η(s) as
∞(cid:88)
n=1
Ω =
Ωn
(2.226)
where
(cid:90) s
(cid:90) s
0
0
Ω1(s) = −
1
2
Ω2(s) =
...
ds1η(s1)
(cid:90) s1
0
ds1
ds2
(cid:2)η(s1), η(s2)(cid:3).
(2.227)
(2.228)
(2.229)
Here, the complications of the time-ordered exponential are moved inside the exponential.
The advantage of this is that truncating Ω at any order will still be anti-Hermitian, and
thus result in a unitary transformation. This is unlike truncating (2.222) which is not
guaranteed to be unitary if any truncations are made. Let’s quickly check that truncating Ω
will still be anti-Hermitian. First, check that η is anti-Hermitian. Starting from U (s)
dU†(s)
ds =
− dU (s)
ds U†(s) with the fact that (AB)† = B†A†
(cid:16)
(cid:17)† =
η†(s) =
U (s)
dU†(s)
ds
dU (s)
ds
U†(s) = −η(s)
(2.230)
which shows that η is indeed anti-Hermitian. To check that the commutators in a term like
75
Ω2 are anti-Hermitian, let’s call A =(cid:2)η(s1), η(s2)(cid:3). Then
A† =(cid:2)η(s1), η(s2)(cid:3)† =
(cid:16)
(cid:17)†
η(s1)η(s2) − η(s2)η(s1)
= η(s2)†η(s1)† − η(s1)†η(s2)†
= (−1)2η(s2)η(s1) − (−1)2η(s1)η(s2)
= −
= −A.
(cid:0)η(s1)η(s2) − η(s2)η(s1)(cid:1)
This proves that the commutator of any two anti-Hermitian operators is itself anti-Hermitian.
Therefore, every term of the Magnus expansion is anti-Hermitian, so truncating at any level
ensures that Ω is anti-Hermitian.
To demonstrate the SRG, let’s consider a small two-level system, represented by the
(2.231)
(2.232)
(2.233)
.
1
1
1 −1
initial Hamiltonian
H = T + V =
76
Let’s try to diagonalize H using the Wegner generator η(s) = [T, H(s)],
1
1
1
0
1 −1
0 −1
η(0) = T H − HT
1
1 −1
1
−
1 −1
−
1
1
1
2
1
1
0
=
=
=
0
0 −1
−1 1
−2 0
= 2iσ2,
(2.234)
and by deﬁnition, Ω(0) = 0. Looking at the recursively deﬁned derivative of Ω
∞(cid:88)
k=0
dΩ
ds
=
Bk
k!
adk
Ω(η),
ad0
adk
Ω(η) = η,
Ω(η) =(cid:2)Ω, adk−1
Ω (η)(cid:3).
(2.235)
At s = 0 we have dΩ
ds |s=0 = η(0), since ad1
Ω = [0, η(0)] = 0. In general, the next step is to
calculate Ω(s) by integrating Eqn. (2.235), and then ﬁnd the transformed Hamiltonian as
H(s) = eΩ(s)H(0)e−Ω(s)
(2.236)
by using the Baker-Campbell-Hausdorf expansion. However, in this simple model, we can just
take the exponential of the Pauli matrices rather than doing a truncated BCH expansion.
77
With this example η and Ω truncate after one term, and will always be antisymmetric
matrices, that is
η(s) = igη(s)σ2
Ω(s) = igΩ(s)σ2.
(2.237)
(2.238)
With this form for Ω, we can look at the exact BCH expansion for H(0) = σ1 + σ3, using
the matrix exponential of a Pauli matrix. In the case where (cid:126)a = aˆn, we have
eia(ˆn·(cid:126)σ) = 1cos(a) + i(ˆn · (cid:126)σ)sin(a).
(2.239)
In our case, we want eΩ = eigσ2 so to get this, ia(ˆn · (cid:126)σ) = igσ2, therefore a = g and ˆn = ˆy.
This leads to
eigσ2 = 1cos(g) + iσ2sin(g)
(2.240)
and for e−Ω just take g → −g. Then use cos(−g) = cos(g) and sin(−g) = sin(g). Thus the
78
transformed Hamiltonian is
H(s) = eigσ2(σ1 + σ3)e−igσ2
=(cid:2)1cos(g) + iσ2sin(g)(cid:3)(cid:2)σ1 + σ3)(1cos(g) − iσ2sin(g)(cid:3)
=(cid:2)1cos(g) + iσ2sin(g)(cid:3)(cid:2)σ1cos(g) − iσ1σ2sin(g) + σ3cos(g) − iσ3σ2sin(g)(cid:3)
=(cid:2)1cos(g) + iσ2sin(g)(cid:3)(cid:2)σ1cos(g) − i2σ3sin(g) + σ3cos(g) + i2σ1sin(g)(cid:3)
=(cid:2)1cos(g) + iσ2sin(g)(cid:3)(cid:2)σ1(cos(g) − sin(g)) + σ3(cos(g) + sin(g))(cid:3)
= σ1(cos(g) − sin(g))cos(g) + σ3(cos(g) + sin(g))cos(g)
+ iσ2σ1(cos(g) − sin(g))sin(g) + iσ2σ3(cos(g) + sin(g))sin(g)
= σ1(cos(g) − sin(g))cos(g) + σ3(cos(g) + sin(g))cos(g)
− i2σ3(cos(g) − sin(g))sin(g) + i2σ1(cos(g) + sin(g))sin(g)
= σ1
(cid:2)cos2(g) − sin2(g) − 2cos(g)sin(g)(cid:3)
(cid:2)cos2(g) − sin2(g) + 2cos(g)sin(g)(cid:3)
(cid:2)cos(2g) − sin(2g)(cid:3) + σ3
(cid:2)cos(2g) + sin(2g)(cid:3),
+ σ3
= σ1
where we used σaσb = δab1 + iεabcσc as well as the trig identity cos2(g) − sin2(g) = cos(2g)
and sin(2g) = 2sin(g)cos(g). This can be generalized slightly for a Hamiltonian of the form
79
H(0) = dσ3 + vσ1, where the full transformed Hamiltonian is
H(s) = eigσ2(vσ1 + dσ3)e−igσ2
=(cid:2)1cos(g) + iσ2sin(g)(cid:3)(vσ1 + dσ3)(cid:2)1cos(g) − iσ2sin(g)(cid:3)
=(cid:2)1cos(g) + iσ2sin(g)(cid:3)(cid:2)vσ1cos(g) − ivσ1σ2sin(g) + dσ3cos(g) − idσ3σ2sin(g)(cid:3)
=(cid:2)1cos(g) + iσ2sin(g)(cid:3)(cid:2)vσ1cos(g) − i2vσ3sin(g) + dσ3cos(g) + i2dσ1sin(g)(cid:3)
=(cid:2)1cos(g) + iσ2sin(g)(cid:3)(cid:2)σ1(vcos(g) − dsin(g)) + σ3(dcos(g) + vsin(g))(cid:3)
= σ1(v ∗ cos(g) − d ∗ sin(g))cos(g) + σ3(d ∗ cos(g) + v ∗ sin(g))cos(g)
+ iσ2σ1(v ∗ cos(g) − d ∗ sin(g))sin(g) + iσ2σ3(d ∗ cos(g) + v ∗ sin(g))sin(g)
= σ1(v ∗ cos(g) − d ∗ sin(g))cos(g) + σ3(d ∗ cos(g) + v ∗ sin(g))cos(g)
− i2σ3(v ∗ cos(g) − d ∗ sin(g))sin(g) + i2σ1(d ∗ cos(g) + v ∗ sin(g))sin(g)
= σ1
(cid:2)v ∗ cos2(g) − v ∗ sin2(g) − 2d ∗ cos(g)sin(g)(cid:3)
(cid:2)d ∗ cos2(g) − d ∗ sin2(g) + 2v ∗ cos(g)sin(g)(cid:3)
(cid:2)v ∗ cos(2g) − d ∗ sin(2g)(cid:3) + σ3
(cid:2)d ∗ cos(2g) + v ∗ sin(2g)(cid:3).
+ σ3
= σ1
(2.241)
With these exact results derived, we can compare the direct SRG integration against
the Magnus formulation. This is shown in ﬁgure 2.1, and it is clear how much numerical
stability is gained from ensuring the unitarity each step with the Magnus expansion. The
Magnus formulation relatively quickly reaches the “machine precision” of the ﬁnite precision
ﬂoating point variables in the calculation, while the direct integration method has error that
is highly dependent on the “time step” δs of the integration.
This 2×2 toy example hides two sources of error that would exist in an actual calculation.
This ﬁrst is that the expressions truncate naturally for η and Ω, which will not happen in
80
Figure 2.1: SRG with direct integration and with the Magnus expansion.
a many-body physics problem, and the second source of error is that the BCH expansion
is done to inﬁnite order in Eqn. (2.241). In a realistic many-body calculation, the terms
in these series are observed to decrease monotonically, and a cutoﬀ tolerance can be used
to minimize the numerical error from truncation. To gain insight into the error generated
from the BCH expansion truncation, the nested commutators can be compared against the
exact expression derived in Eqn. (2.241). In ﬁgure 2.2 this source of truncation error goes
back down to machine precision after about 20 terms, which is pretty substantial, although
machine precision is rarely necessary.
Rather than setting a ﬁxed number of terms in the expansion, it is more typical to specify
an error tolerance so that the BCH expansion can be assuredly not the primary source of
error. This is what is done in ﬁgure 2.3, where the tolerance here is deﬁned to be greater
than the row 0 column 0 element of the next term in the BCH expansion.
By setting = 10−8, we are asserting that we want the smallest eigenvalue (the (0,0)
element of the matrix) to be changing by not more than 10−8 for the next term in the BCH
81
Figure 2.2: Magnus SRG with the exact unitary transformation and with a BCH expansion
truncated after a ﬁxed number of terms
series. There are other tolerances that can be chosen, but looking at the plot, it seems to
work well. As a note, the three tolerances of (10−8, 10−12, 10−16) used (12, 15, 18) terms in
their expansions.
The Magnus expansion has shown to be a great tool for these calculations, but potentially
the greatest gain from this is that any operator can also be computed after the ﬂow. Once
the ﬂow has ﬁnished, it only costs an exact BCH to compute O(s) = eΩOeΩ† along with
H(s) = eΩHeΩ†.
While the 2x2 example was not an exercise in many-body physics, IM-SRG and the use of
the Magnus expansion have seen great success in the calculation of many interesting nuclear
systems [42, 43, 44, 45, 46, 47, 48, 3]. The many-body physics results of this work will focus
on the use of coupled cluster theory which is detailed in Chapter 4.
82
Figure 2.3: Magnus SRG with the exact unitary transformation and with a BCH expansion
truncated after a ﬁxed tolerance is met
83
Chapter 3
Physical Systems
With the mathematical framework of quantum many-body theory as a foundation any
quantum system can be investigated numerically. The many-body Schr¨odinger equation has
proven to be an excellent model for studying nearly any physical system for which the parti-
cles are traveling suﬃciently slower than the speed of light. A wide range of ﬁelds including
atomic physics, quantum chemistry, materials science, and nuclear physics greatly bene-
ﬁt from these theoretical tools, which make studying the mathematical and computational
methods surrounding many-body physics a worthwhile endeavor per se. While the physi-
cal systems introduced in this chapter have applications in answering real world questions,
much of the interest in these systems is theoretical. The pairing model is a simple quantum
system which can be studied analytically and exactly. It is therefore an excellent testing
ground for properties of various many-body methods, and as a system to validate numerical
implementations.
Inﬁnite fermionic matter is important for studying valence electrons in
metals [23], and also for studying the volumetric bulk of neutrons thought to constitute the
crust of neutron stars, or as a model for dense nuclear matter [2, 1, 22].
84
3.1 Pairing Model
The pairing model Hamiltonian ˆH = ˆH0 + ˆV is deﬁned as
(cid:88)
(p − 1)a†pσapσ
pσ
ˆH0 = δ
(cid:88)
pq
ˆV = −
1
2
g
a†p+a†p−, aq−aq+
(3.1)
(3.2)
which represents a quantum system with p levels, each having a spin degeneracy of two.
A common choice for single-particle states are eigenstates of the Hartree-Fock operator,
(ˆu + ˆuHF)|p(cid:105) = p |p(cid:105). In the pairing model, this condition is already fulﬁlled. We deﬁne
the states below the Fermi level as holes and redeﬁne the single-particle energies,
(cid:88)
i
q = hqq +
(cid:104)qi|ˆv|qi(cid:105) .
(3.3)
To be more speciﬁc, let us look at the pairing model with four particles and eight single-
particle states. These states (with δ = 1.0) could be labeled as shown in Table 3.1. The
Hamiltonian matrix for this four-particle problem with no broken pairs is deﬁned by six
possible Slater determinants, one representing the ground state and zero-particle-zero-hole
excitations 0p− 0h, four representing various 2p− 2h excitations and ﬁnally one representing
a 4p − 4h excitation. Ignoring Slater determinants with broken pairs, this problem is then
85
Table 3.1: Single-particle states and their quantum numbers and their energies from
Eq. (3.3). The degeneracy for every quantum number p is equal to two due to the two
possible spin values.
State Label p 2sz E
0
1
2
3
4
5
6
7
1
-1
1
-1
1
-1
1
-1
type
hole
-g/2
hole
-g/2
1-g/2 hole
1-g/2 hole
2
2
3
3
1
1
2
2
3
3
4
4
particle
particle
particle
particle
represented by the Hamiltonian matrix
H =
0
4δ − g −g/2 −g/2
2δ − g −g/2 −g/2 −g/2 −g/2
−0
−g/2
−g/2
−g/2 −g/2
−g/2 −g/2
6δ − g −g/2 −g/2
−g/2 −g/2
−g/2
8δ − g −g/2
10δ − g
−g/2 −g/2 −g/2 −g/2
−g/2 −g/2
6δ − g
0
0
0
0
.
(3.4)
Here, the exact eigenvalues can be found by diagonalizing this small matrix. Additionally, it
is easy to calculate low orders of many-body perturbation theory analytically. This is a very
useful check of the numerical implementation since this analytical expression can also be
used to check our coupled cluster implementation as described in Chapter 4. As a reminder,
the expression for the correlation energy for MBPT(2) is
∆EM BP T 2 =
1
4
(cid:88)
abij
(cid:88)
a** j) and (a > b), however it is often convenient to
just compute these for all i, j, a, b for reasons that will be explained later. For now, let’s just
get a handle on the scaling. If these equations are solved for i, j, a, b that already brings a
computation complexity of O(N 2
number of occupied single-particle states (holes) relative to the Fermi energy, and Np is the
p ) to just loop through the entries of tab
ij , where Nh is the
hN 2
127
number of unoccupied single-particle states. This “big O” notation is used to get a handle on
the scaling of the most expensive term, without worrying about constants of multiplication
or non-leading terms. Within each t-amplitude equation there are some heavy sums, like
(cid:88)
lmde
1
4
de tde
vlm
ij tab
lm,
(4.64)
which brings an additional O(N 2
it is common to ignore the diﬀerence between the number of hole states and particle states
p ) leading to a cost of O(N 4
p ). With the big O notation,
hN 4
hN 2
and to just use N ≈ Nh ≈ Np, and to say that CCD written in this form is an O(N 8) theory
for calculating the ground state energy.
For memory requirements, the primary objects that need to be stored are tab
ij and vab
cd,
which require N 2
hN 2
p and N 4
p number of elements to be stored. In many realistic calculations,
Np ≈ 10 ∗ Nh or even Np ≈ 100 ∗ Nh, meaning that vab
by far the largest object that needs to be stored in this theory.
cd with four particle state indices is
If there are three-body
forces in the calculation, the object wabc
def needs to be stored as well, placing some serious
memory requirements on the calculation. A discussion of the memory requirements and data
structures to handle the requirements of CC are explained in detail in Chapter 5.
It turns out that by reorganizing some terms, we can reduce the computational complexity
of the CCD equations. For example, if we deﬁne the intermediate term
Xlm
ij =
vlm
de tde
ij
(4.65)
(cid:88)
de
128
the term in Eqn. (4.64) can be rewritten with the intermediate term as
(cid:88)
lm
1
4
Xlm
ij tab
lm =
1
4
(cid:88)
lmde
vlm
de tde
ij tab
lm.
(4.66)
By computing the sum over the de indices ﬁrst, storing the intermediate result, then after-
wards computing the sum over the lm indices, we have gone from a scaling of O(N 2
hN 2
p ) to
O(N 2
h + N 2
p ), which is a huge advantage! The only downside to this is that extra storage
must be used to store X, but this is typically small compared to the other storage require-
ments. This means that using intermediates, CCD is actually only a O(N 6) theory, which
is extemely cheap for the accuracy it brings. A systematic way of generating these interme-
diates is to follow the development of diagrams for the coupled cluster eﬀective Hamiltonian
H ≡ ( ˆHN e
ˆT )C as outlined by Shavitt and Bartlett [16].
Any term which is quadratic in t can be done in two steps with an intermediate. However,
by examining which terms can be grouped into their own operators in the CC eﬀective Hamil-
tonian, we can occasionally reuse terms for additional eﬃciency. We deﬁne the intermediates
as
Xa
b = f a
b −
1
2
lmd
1
2
(cid:88)
(cid:88)
(cid:88)
(cid:88)
1
2
del
de
1
2
dl
bd tad
vlm
lm,
detde
vil
jl ,
ij
detde
kl ,
v
dbtda
vil
jl ,
Xi
j = f i
j +
X
ij
kl = v
ij
kl +
Xia
jb = via
jb −
129
(4.67)
(4.68)
(4.69)
(4.70)
and then we can rewrite the CCD equation as
0 =vab
ij + ˆP (ab)
Xa
(cid:88)
d
vab
detde
ij +
(cid:88)
de
1
2
−
(cid:88)
l
Xl
itab
lj
lm − ˆP (ab|ij)
d tdb
ij − ˆP (ij)
(cid:88)
Xlm
ij tab
1
2
lm
(4.71)
(cid:88)
ld
Xlb
idtad
lj ,
which is exactly equivalent to the O(N 8) equations, but reduced down to O(N 6) operations
with minimal additional memory requirements.
The computational details for quantum many-body methods is the primary focus of this
work and speciﬁcally Chapter 5. As a method with polynomial scaling coupled cluster theory
is a great method to investigate computational challenges, since improvements in the data
structures and algorithms implementing these equations can greatly expand what can be
calculated with the method. This is in contrast to exact methods like full conﬁguration
interaction, where it will never outrun the factorial scaling as the problem size increases.
Calculations with a large number of particles (Nh ∼ 100), and a very large basis (Np ∼
104−105) are frequently needed for calculations of interesting physical systems, which require
great care in the performance of the CC implementation. Alternatively, for systems without
as extreme basis size requirements, it may be necessary to use triples (CCDT) ˆT ≈ ˆT2 + ˆT3
or three-body forces ˆWN to achieve the accuracy desired for a calculation, but these bring a
heavy cost.
For a calculation of CCDT with three-body forces, one term for the (cid:104)Φabc
ijk| projected
equations is
1
8
(cid:88)
ef g
wabc
ef gt
ef g
ijk ,
(4.72)
which will result in an O(N 9) scaling theory. This is simply too expensive for anything but
130
the smallest systems. Three-body forces are brought up here to bring attention to other
ways which CC theory can grow prohibitively expensive along with increasing basis size N .
This is a pattern in any many-body method, that increasing the accuracy of a calculation
always has a cost, which dictates which physical systems can be studied and which cannot.
Big O notation can show asymptotically how CC theory will grow, but to consider how
expensive a particular calculation will be, it lacks predictive power. Here the multiplicative
constants (c ∗ N 9) and non-leading terms can be quite important. Chapter 5 will show how
the implementation of the CC equations into a code can vary the multiplicative constant by
up to ﬁve orders of magnitude. This swing in cost is nearly impossible to see just from the
CC equations alone, and so great care should be taken when writing a computer program as
it can determine the viability of a many-body method.
131
Chapter 5
Computational Methodology
The previous chapters have stressed how large and unwieldy a many-body calculation can
be. Even for a modest single-particle basis, the factorial growth of the full Slater determinant
basis becomes quickly impossible for even the largest computers in the world. It was shown in
Chapter 4 that coupled cluster (CC) theory generates expressions for the approximate ground
state energy which have polynomial scaling. This allows CC theory to compute properties
of much larger systems by sacriﬁcing some accuracy of the solution. However, even with
eﬃcient implementations of these equations with intermediate diagrams [65], coupled cluster
theory is still computationally expensive and runs into computational limits for all but very
small physical systems. Modern many-body physics necessarily becomes a computationally
challenging ﬁeld just by the very scale of the problems at hand. This chapter will detail
how the same mathematical expressions on paper can take centuries to compute or seconds
to compute depending on the choice of data structures and algorithms implemented in the
code.
5.1 Code Validation
Before these optimizations are implemented, it is useful to ﬁrst implement the many-body
method equations into code in the most direct translation from mathematics as possible.
Optimizing the code to run faster and compute larger basis sets will increase the number
132
of lines of code substantially, increasing the chance of human error. Once an ineﬃcient,
but correct version of the code is ﬁnished, incremental optimizations moving forward can be
compared to the previously validated solution.
5.1.1 Pairing Model
Here, a simple system like the pairing model described in Chapter 3 is an excellent small
system to check the numerical results. In the case of the simple pairing model it is easy
to calculate ∆EM BP T 2 analytically from Eqn. (2.186), where MBPT2 refers to many-body
perturbation theory that was described in Chapter 2. This is a very useful check of our codes
since this analytical expression can also be used to check our ﬁrst CCD iteration. We restate
this expression here but restrict the sums over single-particle states
∆EM BP T 2 =
1
4
(cid:88)
abij
(cid:88)
(cid:104)ij|ˆv|ab(cid:105)(cid:104)ab|ˆv|ij(cid:105)
ab
ij
=
a**