AB INITIO SIMULATIONS OF LIGHT NUCLEAR SYSTEMS USING
EIGENVECTOR CONTINUATION AND AUXILIARY FIELD MONTE

CARLO

By

Dillon Kyle Frame

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Physics — Doctor of Philosophy

2019

PUBLIC ABSTRACT

AB INITIO SIMULATIONS OF LIGHT NUCLEAR SYSTEMS USING EIGENVECTOR

CONTINUATION AND AUXILIARY FIELD MONTE CARLO

By

Dillon Kyle Frame

Since the beginning of the 20th century, the field of nuclear physics has had a lot of major

developments, from Rutherford’s discovery of the nucleus in 1911, to Chadwick’s discovery

of the neutron in 1932, to the development of quantum mechanics in to 1920’s through the

1940’s, and to the development of quantum chromodynamics (QCD) in the 1970’s.

While there has been a large experimental effort to measure the properties of the nuclear

forces, most models that describe these properties are phenomenological. We know that

QCD is the fundamental theory that describes the interactions in a nucleus, but accurately

simulating even a light nucleus in this framework is beyond the capability of modern su-

percomputers. Clearly, a middle-ground is the approach one needs. One kind of ab initio

method used by our research group is called ”chiral effective field theory” (χ-EFT), where

the nuclear forces are broken down into interactions between only protons, neutrons, and a

force-carrying meson, the pion.

In this work, I will describe two techniques that we have used to simulate the properties

of light nuclear systems. The first is called auxiliary field Monte Carlo, in which we re-write

the interactions between particles as an interaction of a single particle with a fluctuating

background field. The second is a new method called eigenvector continuation, in which we

can compute properties in some simple or easy-to-compute configuration, then accurately

extrapolate to configurations that are not able to be computed directly.

ABSTRACT

AB INITIO SIMULATIONS OF LIGHT NUCLEAR SYSTEMS USING EIGENVECTOR

CONTINUATION AND AUXILIARY FIELD MONTE CARLO

By

Dillon Kyle Frame

In this work, we discuss a new method for calculation of extremal eigenvectors and

eigenvalues in systems or regions of parameter space where direct calculation is problematic.

This technique relies on the analytic continuation of the power series expansion for the

eigenvector around a point in the complex plane.

We start this document by introducing the background material relevant to understand

the basics of quantum mechanics and quantum field theories on the lattice, how we perform

our numerical simulations, and how this relates to the nuclear physics we probe. We then

move to the mathematical formalism of the eigenvector continuation.

Here, we present the foundations of the method, which is rooted in analytic function the-

ory and linear algebra. We then discuss how these techniques are implemented numerically

(with a discussion about the computational costs), and the systematic errors associated with

this method.

Finally, we discuss applications of this method to full, quantum many-body systems.

These include neutron matter, the Bose-Hubbard model, the Lipkin model, and the Coulomb

interaction in light nuclei with LO chiral forces. These systems cover two categories of interest

to the field: systems with a substantial sign problem, or systems that exhibit quantum phase

transitions.

To my family, for being supportive, and for encouraging me to give it my all. To my friends
back home, for keeping me from going insane. To my colleagues and classmates, for helping

me when books weren’t enough. And to my academic advisor, Dean Lee, for always

pushing me to be my best, and for the criticisms that helped me learn from my mistakes.

None of this would have been possible without you all. Thank you.

iii

TABLE OF CONTENTS

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

v

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vi

Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

2.1 Elements of Quantum Mechanics

Chapter 2 Background Material . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 The Schr¨odinger Equation . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Operators, States & Measurements
. . . . . . . . . . . . . . . . . . .
2.2 Lattice Field Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Path Integral Formulation . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Grassmann Path Integral . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Transfer Matrix Formalism . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 Calculating Observables
2.3 Nuclear Forces from Effective Field Theory . . . . . . . . . . . . . . . . . . .
2.3.1 Quantum Chromodynamics (QCD) . . . . . . . . . . . . . . . . . . .
2.3.2 Chiral Effective Field Theory . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Fitting the LECs to NN Scattering Data . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 What is Monte Carlo? . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 The Sign Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.3 Euclidean Time Projection Monte Carlo . . . . . . . . . . . . . . . .
2.4.4 Auxiliary Field Monte Carlo . . . . . . . . . . . . . . . . . . . . . . .
2.4.5 Hybrid Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . .
Shuttle Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.6

2.4 Monte Carlo Methods

Chapter 3 Eigenvector Continuation . . . . . . . . . . . . . . . . . . . . . . .
3.1 Mathematical Formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3
Implementation and Computational Cost . . . . . . . . . . . . . . . . . . . .

Chapter 4 Applications to Nuclear Systems

. . . . . . . . . . . . . . . . . .
4.1 Bose-Hubbard Model
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Simulations of Neutron Matter . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Coulomb Interaction in 12C and 16O . . . . . . . . . . . . . . . . . . . . . .
4.4 EC Solutions of the Lipkin Hamiltonian . . . . . . . . . . . . . . . . . . . .

6
7
7
9
12
12
15
17
20
21
24
26
27
29
30
31
32
36
39
42

44
46
53
65

71
71
75
83
98

Chapter 5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . 110

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

iv

LIST OF TABLES

Table 3.1: Tabulated results from gprof

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

67

Table 3.2: Coefficients of best polynomial fits to computational cost of the EC algo-
rithm, as a function of the size of the lattice L. Since there are L3 points,
we expect the leading term to be proportional to L3. We don’t expect other
terms, but present them for completeness. . . . . . . . . . . . . . . . . . .

Table 3.3: Coefficients of best polynomial fits to computational cost of the EC al-
gorithm, as a function of the number of particles in the system Nf . We
expect a small term proportional to L2, from the matrix elements, and a
large linear term, since we keep track of single-particle wavefunctions.
. .

Table 4.1: Summary of the neutron matter results, for 6 and 14 neutrons, compared
to the direct Monte Carlo calculations. For reference, the relevant nuclear
densities are provided. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 4.2: Summary of the Coulomb results in 12C, computed via direct calculation,
and by eigenvector continuation. Dashed entries are where a singular norm
matrix prevented stable calculations. . . . . . . . . . . . . . . . . . . . . .

Table 4.3: Summary of the Coulomb results in 16O, computed via direct calculation,
and by eigenvector continuation. Dashed entries are where a singular norm
matrix prevented stable calculations. . . . . . . . . . . . . . . . . . . . . .

68

70

83

90

93

v

LIST OF FIGURES

Figure 2.1: Energy scales in nuclear physics. At each step, our resolution changes,
and so do the relevant degrees of freedom. This work follows the blue
arrow; describing the jump from QCD to ab initio methods which rely on
EFTs. Figure from G. F. Bertsch et al [6].
. . . . . . . . . . . . . . . . .

Figure 2.2: The running coupling of QCD. Shown is a plot of the renormalized cou-
pling αs as a function of momentum transfer Q. Shown are the numerous
experimental and numerical calculations of αs, at various Q, and the best
fit value at the Z-boson mass MZ is presented [26]
. . . . . . . . . . . .

Figure 2.3: An organization of the diagrams in Chiral EFT. Each row indicates the
order of the diagram, i.e. how many powers of (Q/Λ) the diagram con-
tributes. The solid lines indicate nucleon lines, and the dashed lines indi-
cate pion lines. The columns represent how many nucleons are involved
in the process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.4: Neutron-proton phase shifts and mixing angles, calculated on the lattice,
and shown as a function of the relative momenta prel. Figure courtesy of
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
N. Li [48].

Figure 2.5: Propagation of the many-body state from (τ = 0) to a final state (τ = τf ).
During the propagation, particles can interact via any of the chiral forces
(LO contact and one-pion exchange shown here.)
. . . . . . . . . . . . .

Figure 2.6: Propagation of some initial single-particle state (τ = 0) to a final state
(τ = τf ). Along the path, the state interacts with the background auxil-
iary fields. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.7: A diagram of the Hubbard-Stratonivich transformation. Here, we can
express the two-body interaction in terms of a one-body interaction with
a fluctuating background field. . . . . . . . . . . . . . . . . . . . . . . . .

23

25

27

29

34

35

37

Figure 2.8: Schematic of Shuttle Algorithm. Courtesy of Bingnan Lu (MSU)

. . . .

43

Figure 3.1: Diagram of convergence region for Eqn. (3.2)

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

47

Figure 3.2: Diagram of Extended Convergence Region, Eqn. (3.5)

. . . . . . . . . .

48

Figure 3.3: Extending the Convergence region to include the point c = v . . . . . . .

48

Figure 3.4: Final region of convergence

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

49

vi

Figure 3.5: D is a finite region of the Riemann surface associated with the vector
|ψ1(c)(cid:105). We assume |ψ1(c)(cid:105) is analytic everywhere except for branch point
singularities at z and ¯z. The dashed blue line shows the path traced around
z. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.6: 2nd Eigenvalue of EC norm matrix, as a function of the step size between
values of the couplings used to generate the EC eigenvectors. Here, the
coupling step size is ∆C = 0.001x.
. . . . . . . . . . . . . . . . . . . . .

50

56

Figure 3.7: Overlap of EC wavefunction with exact ground state. . . . . . . . . . . .

57

Figure 3.8: EC estimate for ground state energy vs EC coupling step size

. . . . . .

58

Figure 3.9: Here, we show the condition number of the norm matrix for 3rd order
EC, in the Lipkin model. A condition number ≈ 1 means the matrix can
be inverted with no loss of precision, whereas a large condition number
indicates numerical instability during the inversion procedure.
. . . . . .

59

Figure 3.10: Steps #1 - #3 for the implementation of the online algorithm . . . . . .

62

Figure 3.11: Step #4 for the implementation of the online algorithm . . . . . . . . . .

63

Figure 3.12: Plots of the perturbative Coulomb energy shift in 4He computed with a
jackknife analysis (in blue), versus the eigenvector continuation estimate
for the Coulomb energy shift (orange). On the top plot, the EC estimate
was done with a jackknife analysis, and on the bottom, it was done using
the online algorithm.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.13: gprof breakdown of the computational time spent in major subroutines
in the MCLEFT code. At 1st order, the overall time spent in the EC
subroutines is about 32%, and at 2nd order, that increases to 56% . . . .

Figure 3.14: Plot of the computational cost/time used by the EC method, as a function
of the length of the edges the lattice L. This was done for a fixed number of
particles (N=4), and a fixed number of timesteps (Ltinner = 20, Ltouter
= 20).
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.15: Plot of the computational cost/time used by the EC method, as a function
of the number of particles in the system. This was done for fixed lattice
volume (L=4), and a fixed number of timesteps (Ltinner = 20, Ltouter =
20).
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.1: Ground state energy for the Bose-Hubbard model, computed by exact
. . . . . . . . . . . . . . . .

diagonalization and by perturbation theory.

64

66

68

69

72

vii

Figure 4.2: Ground state energy of the Bose-Hubbard model, done with up to 5th
. . . . . . . . . . . .

order eigenvector continuation (5 training vectors.)

Figure 4.3: Computation of the lowest two energies in the Bose-Hubbard model. Here,
the lowest two eigenvectors at each training coupling were kept and used
to construct the EC subspace. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.4: LO chiral forces used in neutron matter simulations. Included interactions
include a contact term and one-pion exchange . . . . . . . . . . . . . . .

Figure 4.5: Direct calculation of the ground state energy for six neutrons, shown as a
function of Euclidean time. The red, open squares are the lattice results,
and the solid red line is the best fit showing the exponential decay of the
signal. The dashed red line is the 1σ error band for the fit. . . . . . . . .

Figure 4.6: Direct calculation of the ground state energy for fourteen neutrons, shown
as a function of Euclidean time. The red, open squares are the lattice
results, and the solid red line is the best fit showing the exponential decay
. . . .
of the signal. The dashed red line is the 1σ error band for the fit.

Figure 4.7: Eigenvector continuation results for six neutrons, including up to three
vectors in the training subspace. Open squares show the EC estimates for
the ground state energy for a particular set of vectors used to construct
the subspace. The lines show the best fits of these results.
. . . . . . . .

Figure 4.8: Eigenvector continuation results for fourteen neutrons, including up to
three vectors in the training subspace. Open squares show the EC esti-
mates for the ground state energy for a particular set of vectors used to
construct the subspace. The lines show the best fits of these results. . . .

Figure 4.9: Values of the Coulomb potential on the nz=0 plane. The value at the
origin is chosen to be V (0) = 2αem. This is large enough to give good
results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.10: Direct calculation of Coulomb energy shifts for 12C, at C(cid:12) = 4αem. The
black dotted line shows the perturbative results, and the orange squares
show the results of the direct calculation. We see that for large projection
times the sign problem begins to show up, before we have reached a proper
plateau.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.11: EC Calculation of Coulomb energy shifts in 12C. Again, we show the
perturbative results as a reference. Early on, the EC norm matrix is very
singular, so data is not present before Lt = 40. We note that the 2nd
and 3rd order results are computed at the same projection times; they are
only separated on the plot for visual clarity.
. . . . . . . . . . . . . . . .

73

74

77

79

79

81

82

85

86

88

viii

Figure 4.12: Combined plot showing the direct results and the EC estimates of the
Coulomb energy shifts in 12C. We can immediately see the improvement
gained by the EC method, and the overall resilience to the sign problem.
We note that the three points in each pack are computed at the same
.
projection times; they are only separated on the plot for visual clarity.

Figure 4.13: Computation of Coulomb energy shifts in 12C for arbitrary target cou-
pling. C(cid:12) was varied in the range 0.0αem ≤ C(cid:12) ≤ 10.0αem. Shown in
the bottom figure is a zoom of larger coupling region, showing a clear
deviation from the perturbative results. . . . . . . . . . . . . . . . . . . .

Figure 4.14: Plot of the EC estimate for the ground state energy of 12C, as a function
of the MC step. Shown are also the 1σ error bars at each step. These
calculations were done at a target coupling C(cid:12) = 3αem. We see that after
7000 steps, the results are more-or-less converged, though the 3rd order
. . . . . . . . . . . . . . .
results experience some numerical instability.

Figure 4.15: Direct calculation of Coulomb energy shifts for 16O . . . . . . . . . . . .
Figure 4.16: EC Calculation of Coulomb energy shifts in 16O. For Lt = 10 − 30, the
3rd order results were too unstable to extract reliable data. . . . . . . . .

Figure 4.17: Combined plot showing the direct results and the EC estimates of the
Coulomb energy shifts in 16O. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4.18: Computation of Coulomb energy shifts in 16O for arbitrary target cou-
pling. C(cid:12) was varied in the range 0.0αem ≤ C(cid:12) ≤ 10.0αem. Shown in
the bottom figure is a zoom of larger coupling region. . . . . . . . . . . .

Figure 4.19: Plot of the EC estimate for the ground state energy of 16O, as a function
of the MC step. Shown are also the 1σ error bars at each step. These
calculations were done at a target coupling C(cid:12) = 3αem.
. . . . . . . . .

Figure 4.20: Basic schematic of the Lipkin model. There are N particles distributed
among two angular momentum levels, m = ±1/2, separated by some
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
energy gap 
.

89

91

92

93

94

95

96

97

99

Figure 4.21: Results of the N=14 particle system using perturbation theory, up to 5th
order. It is clear that perturbation theory diverges quickly for couplings
far from 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

Figure 4.22: Results of the N=30 particle system using perturbation theory, up to 5th

order. Convergence here is worse than the 14 particle system.

. . . . . . 100

ix

Figure 4.23: EC Estimates of the five lowest energy levels of the N=14 Lipkin model.
The exact energies, obtained via diagonalizing H, are shown on the top
left. At first order and second order, the estimates are not great, but the
third order estimate (and beyond) are extremely good.

. . . . . . . . . . 101

Figure 4.24: Error between the EC estimates and the exact eigenvalues for N=14 Lipkin
model, shown on a log scale. The four lowest eigenvalues are shown, each
calculated up to 5th order. 10−14 is about the machine level of precision,
so better errors could be obtained by increasing precision.

. . . . . . . . 102

Figure 4.25: EC Estimates of the five lowest energy levels of the N=30 Lipkin model.
The exact energies, obtained via diagonalizing H, are shown on the top
left. The results converge much more slowly, and the improvements at
higher orders are noticeable. . . . . . . . . . . . . . . . . . . . . . . . . . 102

Figure 4.26: Error between the EC estimates and the exact eigenvalues for N=30 Lipkin
model, shown on a log scale. The four lowest eigenvalues are shown,
each calculated up to 5th order. The order-by-order convergence is more
pronounced with this larger number of particles. . . . . . . . . . . . . . . 103

Figure 4.27: EC Estimates of the five lowest energy levels of the N=30 Lipkin model,
with W = 0.4. The exact energies, obtained via diagonalizing H, are
shown on the top left. The 1st, 3rd, and 5th order results are shown,
because convergence is overall slower than the results from W = 0 . . . . 103

Figure 4.28: Error between the EC estimates and the exact eigenvalues for N=30 Lipkin
model, with W = 0.4, shown on a log scale. The four lowest eigenvalues
are shown, each calculated up to 5th order. While the residuals are larger,
the order-by-order convergence is still present. . . . . . . . . . . . . . . . 104

Figure 4.29: Plotted here are the first 20 energy levels for the N=1000 Lipkin Model.
Due to the overall energy scale, we have subtracted out the lowest energy
level. We note that the energy levels pinch together at around C = 1.0,
suggesting a nearby branch point. . . . . . . . . . . . . . . . . . . . . . . 106

Figure 4.30: EC estimates for the lowest 5 eigenvalues of the N=1000 Lipkin model.
We show the energy with the lowest ground state eigenvalue subtracted
out, to show the kink present in the spectrum. The EC results diverge
rapidly past the branch point at C = 1. . . . . . . . . . . . . . . . . . . . 107

Figure 4.31: EC estimates for the lowest 5 eigenvalues of the N=1000 Lipkin model.
We show the energy with the lowest ground state eigenvalue subtracted
out, to show the kink present in the spectrum. The EC results diverge
rapidly past the branch point at C = 1. . . . . . . . . . . . . . . . . . . . 107

x

Figure 4.32: Residuals for the first 4 eigenvalues for the N=1000 Lipkin model. The

first 15 EC orders are shown for each eigenvalue.

. . . . . . . . . . . . . 108

Figure 4.33: EC estimates for the lowest 5 eigenvalues of the N=1000 Lipkin model,
with non-zero W. We show the energy with the lowest ground state eigen-
value subtracted out, to show the kink present in the spectrum. The EC
results diverge rapidly past the branch point at C = 3.

. . . . . . . . . . 108

Figure 4.34: EC estimates for the lowest 5 eigenvalues of the N=1000 Lipkin model,
with non-zero W. We show the energy with the lowest ground state eigen-
value subtracted out, to show the kink present in the spectrum. The EC
results diverge rapidly past the branch point at C = 3.

. . . . . . . . . . 109

Figure 4.35: Residuals for the first 4 eigenvalues for the N=1000 Lipkin model, with

non-zero W. The first 15 EC orders are shown for each eigenvalue. . . . . 109

xi

Chapter 1

Introduction

Eigenvector continuation [28] is a new tool that allows calculations of the properties of quan-

tum many-body systems that are otherwise inaccessible; either through the Monte Carlo sign

problem, or by sheer cost in computation. This technique relies on the analytic continuation

of the power series expansion for the eigenvector around a point in the complex plane. We

assume we have a Hamiltonian that can be written as

H = H0 + cH1

(1.1)

where H0 is some baseline, and H1 is the interacting term, with tunable, complex coupling c

that we wish to compute matrix elements of. In describing the mathematics and numerical

techniques associated with this method, we will be assuming that H is a Hermitian matrix,

and that its eigenvalues and eigenvectors can be accurately calculated for at least a few

values of c.

This work is organized into three main sections. In the first chapter, we begin by in-

troducing the background material relevant to understand the basics of quantum mechanics

and quantum field theories on the lattice, how we perform our numerical simulations, and

how this can be tied back to the nuclear physics we are interested in probing. The most

fundamental QFT, quantum chromodynamics, is the most reasonable place to start this

1

discussion. The QCD Lagrangian [26] is given by

(cid:16)

(cid:88)

q

L =

¯ψq,a

iγµ∂µδab − gsγµtC

µ − mqδab

abAC

(cid:17)

ψq,b − 1
4

µνF Aµν
F A

(1.2)

and it describes an SU(3) color gauge field interacting with quarks and gluons; the building

blocks of nucleons and other fundamental particles. We then move down in energy scale to

something more suited for nuclear physics; chiral EFT. Chiral EFT [56][65] is an effective

field theory built around the chiral symmetry of QCD, which is a symmetry that arises from

the limit of zero quark masses. The mass of the up and down quarks are quite small, so this

limit can be described well by an effective field theory. To do simulations using chiral EFT,

one first needs to decide how many terms in the expansion to keep, and then, one needs to

find a way to constrain the low-energy-constants, or LECs, that each term is responsible for.

We do this by fitting low-energy n-n, n-p, and p-p scattering data, as well as ground state

properties of the deuteron [48][51][52]. All of this input can provide sufficient constraints to

use this framework for heavier nuclear systems.

In Chapter 2, we move to the primary topic of this work; the presentation of the math-

ematical formalism of eigenvector continuation. Here, we present the foundations of the

∞(cid:88)

method, which is rooted in analytic function theory and linear algebra. We take as our
starting point the perturbative expansion of the eigenvectors of H, denoted |ψj(c)(cid:105)

|ψj(c)(cid:105) =

|ψ

(n)
j

(0)(cid:105)cn/n!

(1.3)

n=0

were c is assumed to be in a region in which it converges. Suppose that we are interested in

a value of c outside of this region. Using just perturbation theory, we would be out of luck.

2

However, suppose that we construct a new series about a point |w| < |c|. Then, we could

write

∞(cid:88)
∞(cid:88)

n=0

|ψj(c)(cid:105) =

|ψ

(n)
j

(w)(cid:105) =

|ψ

(n)
j

(w)(cid:105) (c − w)n /n!

|ψ

(n+m)
j

(0)(cid:105)wm/m!

m=0

Combining these, we can obtain

|ψj(c)(cid:105) =

∞(cid:88)

∞(cid:88)

(c − w)n wm|ψ

(n+m)
j

(0)(cid:105)/m!n!

(1.4)

(1.5)

(1.6)

n=0

m=0

This multi-series expansion is the heart of the eigenvector continuation method. In a sense,
we now have an expression for |ψj(c)(cid:105) with a wider area of convergence. Now, from a

numerical standpoint, the analogue of derivatives of the eigenvectors at c = 0 are the finite

differences, i.e. computing the eigenvectors at multiple small values for the couplings.

Once we obtain these eigenvectors, we can form them into a projection operator P and

project the Hamiltonian H at some target coupling C(cid:12) into this lower-dimensional subspace.
We call the projected Hamiltonian M = P†HP , and the norm matrix N = P†P . Then, we

solve the generalized eigenvalue problem

M(cid:126)v = λiN(cid:126)v

(1.7)

The lowest eigenvalue of this system is called the EC estimate for the ground state energy,

and the ground state eigenvector in the full space can be reconstructed from the projection

3

operator

|ψ0(C(cid:12))(cid:105) = P(cid:126)v0

(1.8)

From here, we briefly discuss how these techniques are implemented numerically, with a

discussion about the computational costs.

In Chapter 3, we discuss applications of this method to full, quantum many-body sys-

tems. The first of these is the Bose-Hubbard model, which consists of a system of N bosons

interacting on a 3D lattice. The Hamiltonian of this system is given by

(cid:88)

(cid:104)n,n’(cid:105)

H = −t

a†(n’)a(n) +

U
2

ρ(n) [ρ(n) − 1] − µ

n

n

ρ(n)

(1.9)

(cid:88)

(cid:88)

where a†(n) and a(n) are the creation and annihilation operators, and ρ(n) are the density

operators. This system exhibits a phase transition from a weakly-bound Bose gas to a
strongly-bound cluster around the point U/t = −3.8. This makes perturbative methods

difficult, and a good test for the eigenvector continuation. We also looked at neutron matter.

This difficulty in direct calculations here is that for large projection times, the Monte Carlo

sign problem begins to dominate, so reasonable extraction of the ground state energy is not

possible.

The Lipkin model is another such model were perturbation theory has difficulty. This

model considers a N-fold degenerate two-level system, with N spin-1/2 fermions. The Hamil-

tonian consists of the angular momentum operators Jz,J+, and J−.

H = 
Jz +

V
2

(J 2

+ + J 2−) +

W
2

(J+J− + J−J+)

(1.10)

The states of the system |J, m(cid:105) are described by their total angular momentum J and

4

the projection of that momentum onto the z-axis m. Since no term in the Hamiltonian

mixes different J states, we can focus on a single part of H that has the maximum angular

momentum J = Jmax. For small N V /
, the energy levels of the system are ordered by the

m of that state. For larger N V /
, the energy levels pair up, giving an even and odd parity

state with the same energy.

Finally, we discuss adding the Coulomb interaction into our group’s lattice EFT cal-

culations. This normally is plagued by the sign problem, due to the repulsive nature of

the Coulomb interaction, so this is somewhere that eigenvector continuation would help.

These simulations are done at LO in Chiral EFT, and the Coulomb interaction is added in

non-perturbatively using the Hubbard-Stratonivich transformation

(cid:20)

exp

−Cαt
2

(a†(n)a(n))2

=

(cid:21)

(cid:114) 1

2π

(cid:20)

(cid:90) ∞
−∞ ds exp

−1
2

s(n)2 +

(cid:21)
(cid:112)−Cαts(n)a†(n)a(n)

(1.11)

which allows us to write a two-body interaction in terms of a one-body density interact-

ing with a fluctuating background field. Using this “Coulomb” field, we can simulate the

Coulomb interaction for multiple values of the fine-structure constant αem. Then, the EC

method can give us estimates for the ground state energy at the physical value for αem.

In Chapter 4, we briefly summarize our results, and highlight a few plans that we have

moving forward. The eigenvector continuation technique is proving to be an invaluable tool

for quantum many-body simulations, and we would like to investigate applications to many

other systems and problems.

5

Chapter 2

Background Material

The purpose of this chapter is to provide the information that is essential to the work

presented in this document. Experts in this area should skip ahead to the next chapter for

the new content described in this work.

First, an overview of the essential parts of quantum mechanics are presented. This covers

the discussion of the Schr¨odinger equation and how wavefuctions are described in terms of

bras and kets. Then, we jump into a review of quantum field theory, focusing on the the

lattice methods and effective fields theories used in this work.

Afterwards, we will begin laying down the foundations of lattice field theory. These will

describe a general framework for discretizing the path integrals present in QFT into a form

that can be easily expressed on a computer.

The Monte Carlo methods used by our research group will then be covered in detail,

with accompanying pseduocode when necessary. We use three primary Monte Carlo meth-

ods: Euclidean time projection Monte Carlo, which involves computing powers of a transfer

matrix, auxiliary field Monte Carlo, which allows us to re-write the transfer matrix in terms

of an integral over background fields, and hybrid Monte Carlo, which provides a sophisticated

method for updating these auxiliary fields.

6

2.1 Elements of Quantum Mechanics

When compared to other fields of physics, like classical mechanics or electromagnetism, the

field of quantum mechanics is fairly young, having been developed around the 1920s.

The principles of quantum mechanics can be summarized by two key observations. The

first is that on the smallest scales, charge, light, and energy are all quantized; that is, they

can only come in discrete packets, called “quanta”. This fact was observed during early

experiments like the Millikan oil drop experiment, J.J. Thompson’s e/m experiment, and

measurements of the photoelectric effect.

The second observation was the wave-particle duality exhibited by fundamental particles,

like the electron. This was observed by de Broglie in the double slit experiment, showing

that electrons could constructively and destructively interfere with each other.

We will now begin summarizing the essential information from this field. This will not

even remotely be a complete treatment, so for further reference, see [58], [64], or [68].

2.1.1 The Schr¨odinger Equation

Since it was determined that quantum particles acted like waves, it was a natural assumption

that the equation of motion for these particles was a wave equation. This is what Erwin

Schr¨odinger postulated in 1926.

For light, the wave equation in one dimension is given by

dE 2
dx2 =

dE
dt2

1
c2

(2.1)

with the electric field denoted E(x, t), and with c being the speed of light. This equation has

7

the well-known solution

E(x, t) = E0 cos(kx − ωt)

(2.2)

with k = p/¯h and ω = E/¯h, where p is the momentum of the photon, and E is its energy.

Now, we wish to look at a particle with mass m, and some interaction potential V . From

the de Broglie relations for an electron, we have

¯hω =

¯h2k2
2m

+ V

(2.3)

which suggests the equations of motion will only have one time derivative for the ω, two

space derivatives for the k2 term, and no derivatives for the V term.

The Schr¨odinger equation can then be written as

− ¯h2
2m

∂2Ψ(x, t)

∂2x2 + V (x, t)Ψ(x, t) = i¯h

∂Ψ(x, t)

∂t

(2.4)

The presence of the i in this equation suggests that the solutions to this equation will in gen-

eral be complex-valued functions. These solutions, Ψ(x, t) are referred to as wavefunctions.

It is often easier to factor Ψ(x, t) into two pieces, one that is dependent on x and one

that is dependent on t.

Ψ(x, t) = ψ(x)φ(t)

Substituting this into the Schr¨odinger equation, we get

− ¯h2
2m

φ(t)

d2ψ(x)
dx2 + V (x)ψ(x)φ(t) = i¯hψ(x)

dφ(t)

dt

(2.5)

(2.6)

8

and dividing by ψ(x)φ(t) on both sides, we get

− ¯h2
2m

1

ψ(x)

d2ψ(x)
dx2 + V (x) = i¯h

1

dφ(t)

φ(t)

dt

(2.7)

We notice that each side is a function of only a single variable, so we can treat each side as

a separate differential equation, with a common “separation constant” E, which is the total

energy of the particle. We then get two equations

(2.8)

(2.9)

(cid:34)

(cid:35)

ψ(x) = Eψ(x)

− ¯h2
(cid:20)
2m

(cid:21)

d2
dx2 + V (x)
d
dt

i¯h

φ(t) = Eφ(t)

The first equation is referred to as the time-independent Schr¨odinger equation, often

simplified to Hψ = Eψ, with H being the Hamiltonian of the system. The second equation

will then govern the time evolution of the wavefunction. Just looking at the time-dependent

piece, we see that it is only a 1-D ODE, so it can be solved easily, resulting in

φ(t) = e−iEt/¯h

(2.10)

This is referred to as the time-evolution operator.

2.1.2 Operators, States & Measurements

Now the we have built up a description of the dynamics of a quantum particle in terms

of wavefunctions, it is time to step back and look at how we use these wavefunctions to

calculate properties of the particle.

9

First, we note that the states of the quantum particle live in a complex vector space.

The size and dimension of this space depends on the system being studied; i.e.

for a 1-D

particle in a harmonic oscillator potential can be labelled by n, noting which energy level

the particle is in, and by its x-position.

To describe these vector states, we will use bra-ket notation, developed by Dirac, in which
the state labels are in a ket |α(cid:105) and the complex conjugate is in a bra (cid:104)β|. The bras and

kets are vectors in the full complex vector space, and the labels α and β are the quantum

numbers needed to completely describe a quantum state.

Normal operations of linear algebra can be used when working with these objects: the
inner product of a bra and ket (cid:104)β|α(cid:105) will be a complex number, the norm of a vector is given
by (cid:104)α|α(cid:105), and the outer product is given by |α(cid:105)(cid:104)β|. In the language of linear algebra, the

kets are column vectors, and the bras are row vectors.

If the label on the state is continuous, like say the position x, we still label the states in

terms of bras and kets. In this case, the vector space has an infinite number of dimensions.

We call such a space a Hilbert space.

In general, we assume the different basis states of a given system are orthonormal, mean-

ing they are orthogonal to each other

and normalized

(cid:104)α|β(cid:105) = δαβ

(cid:104)α|α(cid:105) = 1

(2.11)

(2.12)

Since we are working in a linear space, any arbitrary quantum state of a system can be

10

expressed as a linear combination of basis states.

(cid:88)

n

|ψ(cid:105) =

cn|n(cid:105)

This set of basis states |n(cid:105) is called a complete set if

(cid:88)

|n(cid:105)(cid:104)n| = 1

(2.13)

(2.14)

n

The advantage of using this machinery is that we can view operators as matrices that multiply

a ket and give us a new ket.

ˆA|α(cid:105) = ai|α(cid:105)

(2.15)

That is, the whole of quantum mechanics can be turned into solutions of eigenvalue problems.

If we want to perform a measurement, and take the expectation value of some operator

ˆA on some state |α(cid:105), then we compute the quantity

(cid:104)β| ˆA|α(cid:105)

(2.16)

These quantities are usually referred to as “matrix elements”, for the reason that for α, β =
0,··· , N , there will be N × N results of this expectation value.

In our simulations, we consider states with 3 spatial labels, nx, ny, nz, one temporal

label, nt, a label for spin and isospin, S, I, and a label for particle number nf . Then,
our kets will be |nx, ny, nz, nt, S, I, nf >. This is bulky to carry around, so often just gets
simplified to |Ψ(cid:105).

11

2.2 Lattice Field Theory

2.2.1 Path Integral Formulation

One of the fundamental problems in quantum mechanics is the evaluation of how a quantum

state evolves in time. In the Schr¨odinger picture, any quantum state can be written in terms

of energy eigenstates of the Hamiltonian of the system, i.e. solutions of the Schr¨odinger

equation.

ˆH|n(cid:105) = En|n(cid:105)

(cid:88)

|α, t(cid:105) = e−i ˆH(t−t0)|α, t = t0(cid:105) =

e−iEn(t−t0)|n(cid:105)(cid:104)n|α, t = t0(cid:105)

(2.17)

(2.18)

We set t0 = 0, ¯h = c = 1. The wavefunction is this state |α, t > projected into position

n

eigenstates.

ψα (x, t) = (cid:104)x|α, t(cid:105) =

(cid:88)
≡(cid:88)

n

e−iEnt(cid:104)x|n(cid:105)(cid:104)n|α, 0(cid:105)

e−iEntun(x)cn

where un(x) = (cid:104)x|n(cid:105) and

n

cn = (cid:104)n|α, 0(cid:105)

(cid:90)
(cid:90)

=

=

d3x(cid:104)x|n(cid:105)(cid:104)x|α, 0(cid:105)
d3xu∗

n(x)ψα (x, 0)

(2.19)

(2.20)

(2.21)

(2.22)

(2.23)

In quantum field theory, systems are described in terms of fields, φ ((cid:126)x, t), which are

functions of the space-time coordinates. The most fundamental object that describes a

12

system is the path integral, which is defined as

Z =

with Dφ defined as

(cid:90)

(cid:26)

(cid:90) tf

dt S(cid:0)φ, ∂µφ(cid:1)(cid:27)

Dφ exp

i

0

∞(cid:89)

i=0

Dφ =

dφi ((cid:126)x, t)

and S is the action associated with the system

S(cid:0)φ, ∂µφ(cid:1) =

(cid:90)

d3x L(cid:0)φ, ∂µφ(cid:1)

(2.24)

(2.25)

(2.26)

and L is the Lagrangian density (often just called the Lagrangian.) Physically, this path

integral represents a sum over all possible paths that a particle (field) can take, travelling

from one point in time to another, with each path having an associated weight.

It is often easier to work in imaginary time, τ = it, also called Euclidean time. The reason

behind this is that we would like to interpret the weights in Eqn. (2.24) as probabilities,

instead of complex phases. With this substitution, our path integral becomes

(cid:90)

Z =

(cid:26)

−

(cid:90) τf

0

(cid:90)

d3x L(cid:0)φ, ∂µφ(cid:1)(cid:27)

Dφ exp

dτ

(2.27)

Lattice field theory refers to a specific formulation of this path integral, developed by

Wilson [71], in which space-time is discretized on a N+1 dimensional lattice. This is done by

replacing the integrals by sums and by having the fields φµ ((cid:126)x, t) be defined at a discrete set of

points φµ (a(cid:126)n, atnt), where a is the spacing between lattice sites in the N spatial dimensions,

and at is the spacing between lattice sites in the time dimension.

13

The advantage of this method is that it reduces the infinite dimensional path integral

into a finite dimensional, ordinary integral.

(cid:21)

(cid:20)(cid:90)

N(cid:89)

i=0

Z =

dφi

exp

− 1

L3Lt

(cid:88)

(cid:88)

nt

(cid:126)n

L(cid:0)φi ((cid:126)n, nt) , ∂µφi ((cid:126)n, nt)(cid:1) a3at



(2.28)

Thus, using suitable integration methods, this integral can be evaluated, and the dynamics

of the system can be studied.

A free field theory is defined by any involved fields having no interactions; that is, the

lattice action consists of only a kinetic energy term. On the lattice, this corresponds to a

finite difference operation, or a hopping term that allows particles to travel from one site

on the lattice to another. Considering only one spatial dimension, this free lattice action is

defined as:

N(cid:88)

(cid:88)

i=0

x

S = a4

m ˙φ2

i (x) = a4

N(cid:88)

i=0

(cid:88)

x

1
4a

m (φi(x) − φi(x − a))2

(2.29)

An interacting field theory adds potential energy terms between the fields V (|x− x(cid:48)|), which

depends on the distance between lattice sites. In this case, the action has the form

N(cid:88)

(cid:88)

i=0

x

S = a4

m (φi(x) − φi(x − a))2 +

1
4a

a
2

φi(x)V (x − x(cid:48))φi(x(cid:48))

(2.30)

(cid:88)

x,x(cid:48)



If we take the interaction potential as a delta function V (x − x(cid:48)) = Cδ(x − x(cid:48)), then this

reduces to

N(cid:88)

(cid:88)

i=0

x

(cid:20) 1

4a

S = a4

m (φi(x) − φi(x − a))2 +

Ca
2

φ2
i (x)

(cid:21)

(2.31)

14

2.2.2 Grassmann Path Integral

In the previous section, we described the path integral for a scalar field. However, in nuclear

systems, we are concerned with anti-commuting fermionic fields. These are described using

a Grassmann algebra, in which its elements η1, η2, ..., anti-commute

{η1, η1} = 0

(2.32)

These variables have a few interesting properties:

• Linearity

(cid:90)

[af (η) + bg(η)] dη = a

• Integration/Differentiation

(cid:90)

f (η)dη + b

(cid:90)

(cid:21)

f (η)

dη = 0

(cid:90) (cid:20) ∂
(cid:90)
(cid:90)

∂η

1dη = 0

ηdη = 1

g(η)dη

(2.33)

(2.34)

Let ci and c∗
periodic with respect to the spatial lengths of the L3 lattice,

i be anti-commuting Grassmann fields for spin i. The Grassmann fields are

ci((cid:126)n, nt) = ci((cid:126)n + L ˆe1, nt) = ci((cid:126)n + L ˆe2, nt) = ci((cid:126)n + L ˆe3, nt)

(2.35)

15

and anti-periodic along the temporal direction,

ci((cid:126)n, nt) = −ci((cid:126)n, nt + Lt)

The Grassmann path integral is given by

(cid:90)

Z =

DcDc∗ exp [−S(c, c∗)]

(cid:88)

S(c, c∗) = Sfree(c, c∗) + Cαt

ρ↑((cid:126)n, nt)ρ↓((cid:126)n, nt)

(cid:126)n,nt
The action S(c, c∗) consists of the free fermion action

Sfree(c, c∗) =

(cid:88)
(cid:88)

(cid:126)n,nt,i=↑,↓

− h

(cid:126)n,nt,i=↑,↓

l=1,2,3

i ((cid:126)n, nt)ci((cid:126)n, nt + 1) − (1 − 6h)c∗
[c∗
(cid:88)

[c∗
i ((cid:126)n + ˆel, nt)ci((cid:126)n, nt) + c∗

i ((cid:126)n, nt)ci((cid:126)n, nt)]

i ((cid:126)n, nt)ci((cid:126)n − ˆel, nt)]

(2.36)

(2.37)

(2.38)

(2.39)

and an attractive interaction between up and down spins. The Grassmann spin densities ρ↑

and ρ↓ are defined as

ρ↑((cid:126)n, nt) = c∗

↑((cid:126)n, nt)c↑((cid:126)n, nt)

ρ↓((cid:126)n, nt) = c∗

↑((cid:126)n, nt)c↓((cid:126)n, nt)

(2.40)

(2.41)

Now, we wish to extend the definition of the Grassmann path integral to add in a real-

valued, auxiliary field s((cid:126)n, nt). We do this by defining the Grassmann actions

Sj(c, c∗, s) = Sfree(c, c∗) −(cid:88)

Aj [s((cid:126)n, nt)] ·(cid:2)ρ↑((cid:126)n, nt) + ρ↓((cid:126)n, nt)(cid:3)

(2.42)

(cid:126)n,nt

16

and Grassmann path integrals,

(cid:20)(cid:90)

(cid:89)

(cid:21)(cid:90)

djs((cid:126)n, nt)

DcDc∗ exp(cid:2)−Sj(c, c∗, s)(cid:3)

(2.43)

Zj =

(cid:126)n,nt

Here, Aj is a particular functional form involving the auxiliary field. The first possibility

is a Gaussian-integral transformation similar to the original Hubbard-Stratonovich transfor-

mation, denoted by j = 1.

(cid:90)

(cid:90) +∞
−∞ ds((cid:126)n, nt)e

− 1
2 s2((cid:126)n,nt)

d1s((cid:126)n, nt) ≡ 1√
2π

A1 [s((cid:126)n, nt)] =(cid:112)−C1αts((cid:126)n, nt)

(2.44)

(2.45)

(2.46)

(2.47)

Another possibility is a bounded, continuous auxiliary field transformation, denoted by j = 4.

(cid:90) +π

d4s((cid:126)n, nt) ≡ 1
2π

(cid:90)
A4 [s((cid:126)n, nt)] =(cid:112)−C4αt sin (s((cid:126)n, nt))

ds((cid:126)n, nt)

−π

Other transformations are described in [44], but only the two transformations detailed here

are used in our simulations.

2.2.3 Transfer Matrix Formalism

Relating the Grassmann path integral to the transfer matrix, we have

(cid:16)

(cid:17)

Z = Tr

M Lt

(2.48)

17

That is, the path integral can be re-written as the trace of a product of transfer matrix

operators. Let us define the free lattice Hamiltonian as

(cid:88)

(cid:126)n,i=↑,↓

(cid:88)

(cid:88)

(cid:104)

(cid:126)n,i=↑,↓

l=1,2,3

Hfree =

3
m

†
i ((cid:126)n)ai((cid:126)n) − 1

2m

a

†
i ((cid:126)n)ai((cid:126)n + ˆel) + a

†
i ((cid:126)n)ai((cid:126)n − ˆel)

a

and the density operators as

ρ↑((cid:126)n) = a

†
↑((cid:126)n)a↑((cid:126)n)

ρ↓((cid:126)n) = a

†
↓((cid:126)n)a↓((cid:126)n)

(cid:105)

(2.49)

(2.50)

(2.51)

Then, the normal ordered transfer matrix operator is given by

(cid:34)

(cid:88)

(cid:35)

M =: exp

−Hfreeαt − Cαt

ρ↑((cid:126)n)ρ↓((cid:126)n)

:

(2.52)

Like before, we wish to express the Grassmann path integral as a trace of products of

(cid:126)n

transfer matrix operators. Now, the transfer matrix depends on the auxiliary fields,

(cid:20)(cid:90)

(cid:89)

Z =

djs((cid:126)n, nt)

where

M =: exp

(cid:126)n,nt

(cid:40)
−Hfreeαt +

(cid:88)

(cid:126)n

(cid:21)
Tr(cid:2)Mj(s, Lt − 1)··· Mj(s, 0)(cid:3)
Aj [s((cid:126)n, nt)] ·(cid:2)ρ↑((cid:126)n, nt) + ρ↓((cid:126)n, nt)(cid:3)(cid:41)

(2.53)

:

(2.54)

Note that in this transformation, the interaction part of the lattice Hamiltonian in Eqn.

(2.52) involved two-body operators, but in Eq. (2.54), only one-body operators are involved.

Now the the groundwork is laid for generating a baseline interaction using an auxil-

iary field, we discuss a method for adding in extra interactions with a known potential or

18

propagator. We add in a SU(2N)-invariant potential V ((cid:126)n − (cid:126)n(cid:48)) to the lattice Hamiltonian:

Hnew = Hbase +

1
2

: ρi((cid:126)n)V ((cid:126)n − (cid:126)n(cid:48))ρi((cid:126)n(cid:48)) :

(cid:88)

(cid:88)

(cid:126)n,(cid:126)n(cid:48)

i=↑,↓

where

Hbase = Hfree + C

(cid:88)

(cid:126)n

: ρ↑((cid:126)n)ρ↓((cid:126)n) :

(2.55)

(2.56)

as described earlier in this section. We can then define the normal ordered transfer matrix

M as:

M =: exp

−Hfreeαt − Cαt

(cid:88)

(cid:126)n

ρ↑((cid:126)n)ρ↓((cid:126)n) − αt
2

(cid:88)

(cid:126)n,(cid:126)n(cid:48)

(cid:88)

i=↑,↓

 :

ρi((cid:126)n)V ((cid:126)n − (cid:126)n(cid:48))ρi((cid:126)n(cid:48))

(2.57)

We will denote the Fourier Transforms of quantities (or their momentum space versions) in
this section with a tilde. Thus, the Fourier Transform of the potential V ((cid:126)n − (cid:126)n(cid:48)) is given by

˜V (2π(cid:126)ks/L) =

V ((cid:126)n)e2πi(cid:126)n·(cid:126)k/L

(2.58)

(cid:88)

(cid:126)n

(cid:88)

(cid:126)k

(cid:88)

(cid:88)

(cid:126)n,(cid:126)n(cid:48)

nt

We assume that ˜V (2π(cid:126)ks/L) is negative definite. The inverse of the potential, V −1((cid:126)n) is

V −1((cid:126)n) =

1
Ld

e−2πi(cid:126)n·(cid:126)k/L
˜V (2π(cid:126)ks/L)

(2.59)

Now, we can add in an auxiliary field s1, resulting in the lattice action

Snew = Sbase − αt
2

s1((cid:126)n, nt)V −1((cid:126)n − (cid:126)n(cid:48))s1((cid:126)n(cid:48), nt)

(2.60)

where Sbase is given in Eqn. (2.42). The modified transfer matrix, with both auxiliary fields,

19

is then given by

M (s, s1, t) =: exp

(cid:40)
−Hfreeαt +

(cid:88)

(cid:126)n

Aj [s] ·(cid:2)ρ↑ + ρ↓(cid:3) +

Aj(cid:48) [s1] ·(cid:2)ρ↑ + ρ↓(cid:3)(cid:41)

(cid:88)

(cid:126)n

:

(2.61)

where the space and time indices have been suppressed to save line space. For our simulations,

we use j = 1 and j(cid:48) = 4. These Aj [s((cid:126)n, nt)] functionals are the same as described in Section

2.2.2. Substituting in the appropriate Aj[s], our transfer matrix has the form:

sin s((cid:126)n, nt)(cid:2)ρ↑ + ρ↓(cid:3) +

(cid:88)

(cid:112)−Cαt
s1((cid:126)n, nt) ·(cid:2)ρ↑ + ρ↓(cid:3)} :

(cid:126)n

(2.62)

(2.63)

M (s, s1, t) =: exp{ − Hfreeαt +
(cid:88)

(cid:112)−C1αt

(cid:126)n

2.2.4 Calculating Observables

Mentioned earlier, we wish to take as our observable the powers of the transfer matrix. This

Euclidean time projection amplitude ZN,N (t) is defined as

(cid:20)(cid:90)

ZN,N (t) ≡ (cid:89)

(cid:126)n,nt

(cid:21)

djs((cid:126)n, nt)

(cid:104)Ψ0,free

N,N |Mj(s, t)··· Mj(s, 0)|Ψ0,free
N,N (cid:105)

(2.64)

As a result of the normal ordering, the transfer matrix only consists of single-particle oper-

ators acting with the background auxiliary field. There are no direct interactions between

particles. Thus,

(cid:104)Ψ0,free

N,N |Mj(s, t)··· Mj(s, 0)|Ψ0,free

N,N (cid:105) =(cid:2)det Mj(s, t)(cid:3)2

(2.65)

20

where

(cid:2)Mj(s, t)(cid:3)

k,k(cid:48) = (cid:104)(cid:126)p(cid:48)

k|Mj(s, t)··· Mj(s, 0)|(cid:126)pk(cid:105)

(2.66)

for matrix indices k, k(cid:48) = 1,··· , N .

|(cid:126)pk(cid:105) and |(cid:126)p(cid:48)

k(cid:105) re the single-particle momentum states

comprising the initial Slater determinant state. The single particle interactions in the transfer

matrix are the same for up and down spin particles, which results in the determinant being

squared in Eqn.

(2.65). This matrix is real-valued, so the square of the determinant is

non-negative and there is no sign problem.

To extract the ground state energy, we compute ZN,N (t) and ZN,N (t− αt), and take the

ratio:

EN,N (t) =

Λ
αt

log

(cid:18) ZN,N (t − αt)

(cid:19)

ZN,N (t)

For large t, this quantity converges to the ground state energy:

E0

N,N = lim

t→∞ EN,N (t)

(2.67)

(2.68)

2.3 Nuclear Forces from Effective Field Theory

Nuclear physics is field that spans multiple scales; from the level of quarks and gluons, to

the level of neutrons and protons, to the scale of collectivity in large nuclei. This hierarchy

is summarized in Fig. [2.1]. This hierarchy can be broken into energy scales, which represent

the “resolution” of the associated model. At each level, the properties being investigated

don’t tend to depend on the higher energy scales; if you care about the rotational bands

in Lead, it doesn’t really matter what quark #534 is doing. The higher energy degrees of

freedom are “integrated out”, such that they don’t have any influence in the dynamics of

21

the system.

In this section, we will try to summarize the stepping stones we need to get from the

fundamental interactions in quantum chromodynamics to chiral EFT, which is what we use

in our lattice simulations. Apart from some beyond-the-standard-model theories, QCD is

the most fundamental description of nucleons, in terms of their constituent particles, the

quarks and gluons. However, systems of single protons or neutrons are at the cutting-edge

for lattice QCD calculations. It will be a long time before supercomputers and the numerical

algorithms are able to push the limits into nuclear systems.

22

Figure 2.1: Energy scales in nuclear physics. At each step, our resolution changes, and so
do the relevant degrees of freedom. This work follows the blue arrow; describing the jump
from QCD to ab initio methods which rely on EFTs. Figure from G. F. Bertsch et al [6].

23

2.3.1 Quantum Chromodynamics (QCD)

Quantum Chromodynamics (QCD) is an SU(3) gauge field theory that describes the strong

interactions between quarks and gluons.

It was developed throughout the 1970s, as an

application of Yang-Mills theory.

The QCD Lagrangian is given by

(cid:16)

(cid:88)

q

L =

¯ψq,a

iγµ∂µδab − gsγµtC

µ − mqδab

abAC

(cid:17)

ψq,b − 1
4

µνF Aµν
F A

(2.69)

where repeated indices are summed over. The γµ are the Dirac γ-matrices. The ψq,a is the

quark field spinors, for quark flavour q = u, d, s, c, t, orb, with corresponding masses mq, and
a color index a = 1,··· , Nc = 3, i.e. quarks come with three “color charges.” These quark

fields are therefore elements of fundamental representation of SU(3). The AC
fields, where µ is once again the Dirac index, and C = 1,··· , N 2

µ are the gluon
c − 1 = 8 runs over the

8 kinds of gluons. These gluon fields transform under the adjoint representation of SU(3).

The tC

ab are the Gell-Mann matrices, defined by the relation

(cid:104)
tA, tB(cid:105)

= ifABC tC

(2.70)

where fABC are the structure factors of SU (3). Finally, the field strength tensor F A

µν is

defined by

F A
µν = ∂µAA

ν − ∂νAA

µ − gsfABC AB

µ AC
ν

(2.71)

The coupling of the strong interaction is gs. These enter into the three fundamental

diagrams of QCD: a quark-antiquark-gluon vertex, a three-gluon vertex, and a four-gluon

vertex. These diagrams contribute factors of gs, gs, and g2

s , respectively.

24

Figure 2.2: The running coupling of QCD. Shown is a plot of the renormalized coupling αs
as a function of momentum transfer Q. Shown are the numerous experimental and numerical
calculations of αs, at various Q, and the best fit value at the Z-boson mass MZ is presented
[26]

One peculiarity of QCD is that quarks and gluons cannot be observed as free particles;

instead they must be observed as color-neutral (color singlet) states. These can form mesons

(a quark-antiquark) and baryons (three quarks). This property is known as “confinement.”

Another interesting feature of QCD is that of asymptotic freedom; the strength of the

strong interaction gets weaker at higher energies, and stronger at lower energies. This running

coupling is shown in Fig. [2.2]. For low energies, QCD is non-perturbative, so it is often the

case that computations are done at higher values for the coupling, the extrapolated down to

the “physical point”.

25

2.3.2 Chiral Effective Field Theory

Chiral effective field theory (χ-EFT) and chiral perturbation theory (χ-PT) are built from a

low energy effective field theory that was designed to be consistent with the chiral symmetry

of the QCD Lagrangian. This chiral symmetry, simply put, is the limit where the quark

masses go to zero.

Construction of the chiral Lagrangian starts with considering all terms that are consistent

with the approximate chiral symmetry of QCD. This symmetry can be spontaneously broken,

giving rise to the Goldstone bosons, the pions.

In this low-energy EFT, the relevant degrees of freedom move from the quarks and gluons

to just the hadrons. In χ-EFT, these are taken to be the nucleons. Then, one can work in

terms of a power counting scheme in terms of some momentum scale.

If one looks at a given diagram, the behaviour can be analysed by rescaling the external
momentum lines pi → tpi and the quark masses mi → t2mi. Then, Weinberg’s power
counting shows that a given digram scales as tD for D ≥ 2 is the number of vertices in the
diagram. For a given D, this tells us to what order in the expansion of Lχ we need.

These expansions tend to work for momentum smaller than the chiral symmetry breaking
scale ΛχSB ≈ 4πf0, where f0 = 93 MeV is the pion decay constant. Additionally, the other
scale where χ-PT tends to break down is near the ρ mass, which is ≈ 780 MeV.

Shown in Fig. [2.3] is the organisation of diagrams in terms of this power counting scheme.

In the rows are all diagrams which contribute a given power of (Q/Λ). In the columns are

diagrams which consist of N external nucleon lines.

From this, we can see that the three-nucleon force does not make an appearance until

N2LO, which contributes at the (Q/Λ)3 level.

26

Figure 2.3: An organization of the diagrams in Chiral EFT. Each row indicates the order
of the diagram, i.e. how many powers of (Q/Λ) the diagram contributes. The solid lines
indicate nucleon lines, and the dashed lines indicate pion lines. The columns represent how
many nucleons are involved in the process.

For more information on chiral perturbation theory or chiral effective field theory, see

[65] and [56].

2.3.3 Fitting the LECs to NN Scattering Data

Chiral EFT has several low-energy constants, or LECs, that need to be constrained before

any calculations can be done. These are usually done by looking at low-energy scattering

experiments between nucleons, such as n-n, n-p, or p-p scattering, or by binding energies of

light systems, like the deuteron or alpha particle.

The first step in determining the LECs is to express the relevant operators in terms of

angular momentum states. Then, the operators that act on a given partial wave can have

their own particular LEC. By looking at the phase shifts, these can be fitted directly with

lattice calculations. This section follows the discussion in [48]. For a full discussion, please

27

see this reference.

To start determining the phase shifts on the lattice, we need to re-write the full wave-

function in terms of spherical harmonics. This is done by

|r(cid:105)L,Lz =

YL,Lz (ˆr(cid:48))δ(|r(cid:48)| − r)|r(cid:48)(cid:105)

(2.72)

(cid:88)

r(cid:48)

where r is a fixed distance from the origin, and the sum runs over all lattice points. The

delta function counts only the lattice points that are a distance r from the origin.

With this definition of the radial wavefunction, the Hamiltonian can be reduce to a form

that contains only a single variable, r. Then, the phase shifts and mixing angles can be

extracted in the region of vanishing potentials. Here, the radial wavefunction can be written

in terms of spherical Bessel functions

|r(cid:105)L,Lz ≈ ALh−

L (kr) + BLh+

L (kr)

(2.73)

where k =

√
2µE, µ is the reduced mass, and E is the energy. The scattering coefficients

AL and BL are related to the S-matrix by

BL = SLAL

(2.74)

with SL = exp(2iδL), where δL is the phase shift. This phase shift is determined by

(cid:18) BL

(cid:19)

AL

δL =

1
2i

log

(2.75)

Fitting the phase shifts is done using a least-squares minimization. Shown in Fig. [2.4]

28

Figure 2.4: Neutron-proton phase shifts and mixing angles, calculated on the lattice, and
shown as a function of the relative momenta prel. Figure courtesy of N. Li [48].

are the results of fits of the phase shifts and mixing angles for several partial waves, done on
a lattice with lattice spacing a = 1.97 fm, for the energy range Elab ≤ 50 MeV.

2.4 Monte Carlo Methods

In this section, we will discuss the main methods that we use to evaluate these discretized

path integrals.

29

2.4.1 What is Monte Carlo?

Monte Carlo algorithms are a type of algorithm that involves repeatedly generating random

numbers, and using them in some meaningful way. They were named after the city of Monte

Carlo in Monaco, since these algorithms have the gambling/rolling dice feel to them.

One of the primary uses of Monte Carlo is in numerically integrating a function. If we

consider a function of one variable f (x), over some closed interval [ab], then we can view the

integral of this function over this interval

(cid:90) b

a

f (x)dx

(2.76)

as the area under the curve. All that we need to do is pick points in the domain and range

of f (x), and check to see if the point is above or below f (x). Then, the integral can be

evaluated as

f (x)dx = lim
N→∞

N

(2.77)

The other primary use of Monte Carlo algorithm is in sampling probability distributions.

For a given continuous probability distribution P(x) and operator O(x)

(cid:90) b

a

(cid:90)

(cid:18) Nbelow

(cid:19)

N(cid:88)

i=1

(cid:104)O(cid:105) =

P(x)O(x)dx = lim
N→∞

1
N

P(xi)O(xi)

(2.78)

That is, any expectation value which involves integrating a quantity over some weight distri-

bution can be expressed as a discrete sum over a finite number of samples pulled from that

distribution.

In lattice field theory, this probability distribution is the Boltzmann distribution e−S ,

and the integration variables are the fermion fields. Now, since there is an infinite sum

30

for each variable being integrated, direct sampling would never work for these lattice field

theories. For this, we need more specialized algorithms.

2.4.2 The Sign Problem

Before we hop into the discussion of the Monte Carlo algorithms we use in our simulations,

we first make a few remarks about the sign problem.

The sign problem/fermion sign problem/signal-to-noise problem is a common issue that

arises in the simulations of fermions, and it is due to the fact that the probability distribution

is not necessarily positive-definite. If we consider the expectation values of some operator

(cid:82) DφA[φ]p[φ]
(cid:82) Dφp[φ]

(cid:104)A(cid:105) =

with bosonic fields φ, and the probability weight defined by

p[φ] = det (M [φ]) exp (−S[φ])

(2.79)

(2.80)

which contains the fermion determinant and the Boltzmann weighting factor. When p[φ] is

positive for all φ, this integral is well-defined, and can be evaluated numerically.

However, if p[φ] is ever non-positive (usually due to a complex S), then the integral

becomes highly oscillatory, can cannot reliably be determined numerically. This can be

shown if we break the probability weights into a modulus and a phase

p[φ] = ρ[φ] exp(iθ[φ])

(2.81)

31

and then we move the phase into the definition of the expectation value.

(cid:104)A[φ] exp(iθ[φ])(cid:105)

(cid:104)exp(iθ[φ])(cid:105)

=

(2.82)

(cid:82) Dφρ[φ] exp(iθ[φ])
(cid:82) Dφρ[φ] exp(iθ[φ])

(cid:104)A(cid:105) =

This ρ[φ] is now strictly real and positive-definite, so can always be interpreted as a proba-

bility weight. This phase factor now encodes all of the complex behaviour of the observables
and fields. For exp(iθ[φ]) ≈ 1, there is no sign problem, and this expectation value is well-
defined. If, however, exp(iθ[φ]) → 0, both the numerator and denominator become small,

leading to wildly fluctuating estimates for the expectation value.

This problem is an NP-Hard problem, meaning that a general solution to it is likely

not possible (see P = NP). Instead, the best we can do is to find ways to circumvent

the problem. There are multiple ways to minimize this problem in the context of lattice

simulations; meron cluster algorithms [11], the complex Langevin method [1], a density of

states method [43], pseudofermions [29], Canonical methods [3, 14], analytic continuation

from imaginary chemical potentials [15], fermion bags [9], re-weighting methods [27], sign-

optimized manifolds [2], or thimble methods [12, 13]. Each has its own particular formulation,

advantages, and disadvantages. In this work, we will show how the eigenvector continuation

method can avoid this problem from becoming too much of a problem.

2.4.3 Euclidean Time Projection Monte Carlo

If we are interested in studying the dynamics of a quantum system on the lattice, then it

would be natural to look at the time evolution of a system. However, in quantum mechanics,

32

the time evolution operator is simply a phase factor

ˆT (ti − tj) = e

−i ˆH(ti−tj )

(2.83)

It is then a “natural” extension to consider the case t → iτ , which is a Wick rotation for the

time variable. This new “time” variable τ is referred to as the Euclidean time, or projection

time. The effect of this transformation is it turns the time projection operator into one that

can be calculated easily, as it is a simple exponential decay

ˆTE(τi − τj) = e

− ˆH(τi−τj )

(2.84)

In Monte Carlo simulations, we keep track of the Euclidean time projection amplitude

ZN,N (τ ) ≡ (cid:89)

(cid:20)(cid:90)

(cid:21)

djs((cid:126)n, nt)

(cid:104)Ψ0,free

N,N |Mj(s, t)··· Mj(s, 0)|Ψ0,free
N,N (cid:105)

(2.85)

(cid:126)n,nt

Combining this method with the auxiliary fields, we can break the propagation involving

many-body interactions (shown in Fig.

[2.5]) into the propagation of single particle states

through interactions with fluctuating background fields. This is shown in Fig. [2.6].

We now define an online energy expectation value that depends on the Euclidean time τ

EN,N (τ ) =

1
αt

log

ZN,N (τ − at)

ZN,N (τ )

Looking at the spectral decomposition of ZN,N , we see

ZN,N (τ ) =

|ck
N,N|2e

−Ek

N,N τ

(cid:88)

k

33

(2.86)

(2.87)

Figure 2.5: Propagation of the many-body state from (τ = 0) to a final state (τ = τf ).
During the propagation, particles can interact via any of the chiral forces (LO contact and
one-pion exchange shown here.)

34

Figure 2.6: Propagation of some initial single-particle state (τ = 0) to a final state (τ = τf ).
Along the path, the state interacts with the background auxiliary fields.

35

For the limit of τ → ∞, the lowest energy eigenstate is the primary contributor to the

amplitude, and all the higher energy states are exponentially suppressed.

EN,N (τ ) ≈ E0

N,N +

(cid:88)

k(cid:54)=0

(cid:32)|ck
N,N|2
N,N|2
|c0

(cid:33) e

 e

(Ek

N,N−E0
N,N )at − 1
αt

−(Ek

N,N−E0

N,N )τ

(2.88)

This can be further simplified by noting that the difference between the energy eigenstates

N,N ) is usually smaller than the energy cutoff scale 1/αt. This can simplify then

(Ek

N,N − E0

to

EN,N (τ ) ≈ E0

N,N +

(cid:33)(cid:16)

(cid:32)|ck
N,N|2
|c0
N,N|2

(cid:88)

k(cid:54)=0

(cid:17)

N,N − E0
Ek

N,N

−(Ek

N,N−E0

N,N )τ

e

With this, we say that the ground state energy E0

N,N is the large τ limit

E0

N,N = lim

τ→∞ EN,N (τ )

(2.89)

(2.90)

2.4.4 Auxiliary Field Monte Carlo

Now that we have shown how to evaluate the ground state energy by propagating single-

particle states through Euclidean time, we now wish to show how the auxiliary fields for the

transfer matrices are generated.

The auxiliary fields are defined used a Gaussian integral transformation, known as the

(cid:20)

Hubbard-Stratonovich [35] transformation.

exp

−Cαt
2

(a†(n)a(n))2

=

(cid:21)

(cid:114) 1

2π

(cid:90) ∞
−∞ ds exp

(cid:20)

−1
2

s(n)2 +

(cid:21)
(cid:112)−Cαts(n)a†(n)a(n)

(2.91)

This transformation allows us to decompose two-body interactions into the integral over

one-body interactions with a normally-distributed, background field. The only price is the

36

Figure 2.7: A diagram of the Hubbard-Stratonivich transformation. Here, we can express
the two-body interaction in terms of a one-body interaction with a fluctuating background
field.

addition of an integral at each lattice site for each interaction to be written in this manner.

However, this is where Monte Carlo is well-suited.

For simple contact interactions, the auxiliary field can be directly sampled from a normal

distribution

P = exp

(cid:19)

(cid:18)

−1
2

s(n)2

However, for more complicated interactions, like Coulomb or one-pion exchange, we have to

work harder. Here, we take as our distribution

P = exp (−S) = exp

s1((cid:126)n, nt)V −1((cid:126)n − (cid:126)n(cid:48))s1((cid:126)n(cid:48), nt)

(2.93)

(2.92)



−(cid:88)

(cid:126)n,(cid:126)n(cid:48)

where this action is the same as in Eqn. (2.60), except we have removed the constant factor

in front as well as the sum over time steps.

This form of the action is difficult to sample by itself, so we factor it into a Gaussian-like

form. Defining a t((cid:126)n, nt) variable, we would like

(cid:88)

(cid:126)n,(cid:126)n(cid:48)

s1((cid:126)n, nt)V −1((cid:126)n − (cid:126)n(cid:48))s1((cid:126)n(cid:48), nt) =

[t((cid:126)n, nt)]2

2

(2.94)

37

so that the probability distribution simplifies to

P = exp (−S) = exp

(cid:32)− [t((cid:126)n, nt)]2

(cid:33)

2

By inspection, we see that if we define t((cid:126)n, nt) as

t((cid:126)n, nt) =

V −1/2((cid:126)n − (cid:126)n(cid:48))s1((cid:126)n(cid:48), nt)

(cid:88)

(cid:126)n(cid:48)

(2.95)

(2.96)

we reproduce the original distribution. If we invert this equation to solve for s1((cid:126)n(cid:48), nt), we

have

(cid:88)

(cid:126)n(cid:48)

s1((cid:126)n, nt) =

V +1/2((cid:126)n − (cid:126)n(cid:48))t((cid:126)n(cid:48), nt)

(2.97)

This equation provides a simple prescription for generating a new s1((cid:126)n, nt). We just need

to generate a set of normally-distributed numbers t((cid:126)n, nt), and plug those in. However, this

is a slow process, so using a convolution can speed this updating up tremendously.

If we define ˜u((cid:126)k) to be,

(cid:16) ˜V +1/2((cid:126)k)

˜u((cid:126)k, nt) =

(cid:17) ∗(cid:16)

(cid:17)

˜t((cid:126)k, nt)

(2.98)

where ∗ indicates a convolution, then the auxiliary field s1((cid:126)n, nt) is simply

s1((cid:126)n, nt) = u((cid:126)n, nt)

(2.99)

where the tilde represents the momentum space representations of each variable. These are

obtained through the Fourier transform. With these formulas, the algorithm for updating

the field is then:

38

Step #1: Generate ˜V +1/2((cid:126)k). This only needs to be done once, at the beginning of the code.

Step #2: For each (cid:126)n, nt, generate t((cid:126)n, nt) according to Eqn. (2.96).

(cid:27)

P [t((cid:126)n, nt)] ∝ exp

−1
2

[t((cid:126)n, nt)]2

Step #3: For each time step nt, generate the quantities

(cid:26)

(cid:88)

(cid:126)n

(cid:88)

(cid:126)k

˜t((cid:126)k, nt) =

1

L3/2

t((cid:126)n, nt)e2πi(cid:126)n·(cid:126)k/L

˜u((cid:126)k, nt) = ˜V +1/2((cid:126)k)˜t((cid:126)k, nt)

s1((cid:126)n, nt) =

1

L3/2

˜u((cid:126)k, nt)e−2πi(cid:126)n·(cid:126)k/L

(2.100)

(2.101)

(2.102)

(2.103)

Once observables are computed using this auxiliary field configuration, go back to step

#2 and generate a new configuration. Since the fields are being directly sampled from

the respective probability distribution, there is no need for any importance sampling or

Metropolis accept/reject stages.

This technique is very straight-forward, but usually not the most optimal way to generate

field configurations. For this, we turn to more sophisticated updating algorithms, like Hybrid

Monte Carlo.

2.4.5 Hybrid Monte Carlo

Recall, our goal is to calculate the ratio

ZN,N (t − at)

ZN,N (t)

39

(2.104)

using a Markov chain Monte Carlo method. For t = Ltat, auxiliary field configurations for

the base interaction are sampled according to the weight

exp(cid:8)−Uj(s) + 2 log(cid:2)| det Mj(s, t)|(cid:3)(cid:9)

where

and

(cid:104)

s((cid:126)n, nt)2(cid:105)

(cid:88)

(cid:126)n,nt

1
2

U1(s) =

U4(s) = 0

(2.105)

(2.106)

(2.107)

This importance sampling is handled by the Hybrid Monte Carlo (HMC) method [17]. This

involves computing molecular dynamics trajectories using the Hamiltonian

Hj(s, p) =

1
2

(cid:88)

(cid:126)n,nt

[p((cid:126)n, nt)]2 + Vj(s)

(2.108)

where p((cid:126)n, nt) is the conjugate momentum of s((cid:126)n, nt) and

Vj(s) = Uj(s) − 2 log(cid:2)| det Mj(s, t)|(cid:3)

The algorithm for updating the field is as follows:

Step #1: Select an arbitrary initial configuration s0.

Step #2: Select a configuration p0 according to a normal distribution.

(cid:104)

P

(cid:105) ∝ exp

p0((cid:126)n, nt)

(cid:26)

(cid:104)

−1
2

p0((cid:126)n, nt)

(cid:105)2(cid:27)

(2.109)

(2.110)

40

Step #3: For each (cid:126)n, nt, let

˜p0((cid:126)n, nt) = p0((cid:126)n, nt) − 
step
2

(cid:20) ∂Vj(s)

(cid:21)

∂s((cid:126)n, nt)

s=s0

for some small positive 
step. The derivative of Vj(s) is computed using

∂Vj(s)
∂s((cid:126)n, nt)

=

=

∂Uj(s)
∂s((cid:126)n, nt)

−

∂Uj(s)
∂s((cid:126)n, nt)

− 2

(cid:88)

k,l

M−1

j

det Mj

2

(cid:88)

(cid:104)

k,l

(cid:3)

∂(cid:2)Mj
(cid:3)

kl
∂s((cid:126)n, nt)

∂ det Mj

(cid:3)
∂(cid:2)Mj
∂(cid:2)Mj
(cid:105)

kl

kl
∂s((cid:126)n, nt)

lk

Step #4: For steps i = 0, 1,··· , Nstep − 1, let

si+1((cid:126)n, nt) = si((cid:126)n, nt) + 
step ˜pi((cid:126)n, nt)
˜pi+1((cid:126)n, nt) = ˜pi((cid:126)n, nt) − 
step

(cid:20) ∂Vj(s)

∂s((cid:126)n, nt)

(cid:21)

s=si+1

for each (cid:126)n, nt.

Step #5: For each (cid:126)n, nt, let

pNstep((cid:126)n, nt) = ˜pNstep((cid:126)n, nt) − 
step
2

(cid:20) ∂Vj(s)

(cid:21)

∂s((cid:126)n, nt)

s=sstep

Step #6: Select a random number r ∈ [0, 1). If,

(cid:104)−H(sNstep, pNstep) + H(s0, p0)

(cid:105)

r < exp

(2.111)

(2.112)

(2.113)

(2.114)

(2.115)

(2.116)

then set s0 = sNstep. Otherwise, leave s0 as is.

Once observables are computed using this auxiliary field configuration, go back to step #2

41

and start a new trajectory. The observable we calculate at each configuration is

(cid:2)det Mj(s, (Lt − 1)αt)(cid:3)2
(cid:2)det Mj(s, Ltαt)(cid:3)2

Oj(s, Ltat) =

(2.117)

by taking the ensemble average of this observable over many HMC trajectories, we obtain

ZN,N (t − αt)

ZN,N (t)

(2.118)

2.4.6 Shuttle Algorithm

Now that we have algorithms for generating auxiliary field configurations and propagating

the states through Euclidean time, we are ready to put everything together. In older codes,

we used to generate all of the auxiliary fields first, then propagate the initial wavefunctions

through Euclidean time, and then calculate the expectation values of observables at the

middle time-step. This order is not very efficient, and wastes a lot of computational time.

A new algorithm we began to use is the Shuttle algorithm, discussed briefly in [51]. This

algorithm instead updates the time steps one at the time, sweeping back and forth from

τ = 0 to τ = Ltat. The schematic of this method is shown in Fig. (2.8).

The global update here is the older algorithm, and the local update is this new algorithm.

On the diagram, the red box shows the ”current” time-step being updated. As this red box

moves left and right, new auxiliary fields are generated, using either HMC for the full transfer

matrix steps, or the simpler Metropolis updates for the SU(4)-invariant time steps.

When the red box reaches the middle time-step, the expectation value of all of the

operators are computed. And finally, when the red box reaches the first or last time step,

the initial states are re-updated, and the sweep direction is reversed. Using this algorithm,

42

Figure 2.8: Schematic of Shuttle Algorithm. Courtesy of Bingnan Lu (MSU)

we can minimize the wasted time during the Monte Carlo simulations.

43

Chapter 3

Eigenvector Continuation

In this chapter, we discuss the main topic of this dissertation; the Eigenvector Continuation

technique. This method begins with some Hermitian Hamiltonian matrix H

H = H0 + cH1

(3.1)

with the eigenvectors denoted |ψj(c)(cid:105). Generating a power series expansion for |ψj(c)(cid:105), we

get

∞(cid:88)

n=0

|ψj(c)(cid:105) =

|ψ

(n)
j

(0)(cid:105)cn/n!

(3.2)

were c is assumed to be in a region in which this series converges. Suppose that we are

interested in a value of c outside of this region. Using just perturbation theory, we would be
out of luck. However, suppose that we construct a new series about a point |w| < |c|. Then,

we could write

∞(cid:88)
∞(cid:88)

n=0

|ψj(c)(cid:105) =

|ψ

(n)
j

(w)(cid:105) =

|ψ

(n)
j

(w)(cid:105) (c − w)n /n!

|ψ

(n+m)
j

(0)(cid:105)wm/m!

m=0

(3.3)

(3.4)

44

Combining these, we can obtain

|ψj(c)(cid:105) =

∞(cid:88)

∞(cid:88)

(c − w)n wm|ψ

(n+m)
j

(0)(cid:105)/m!n!

(3.5)

n=0

m=0

This multi-series expansion is the heart of the eigenvector continuation method. In a sense,
we now have an expression for |ψj(c)(cid:105) with a wider area of convergence. Now, from a

numerical standpoint, the analogue of derivatives of the eigenvectors at c = 0 are the finite

differences, i.e. computing the eigenvectors at multiple small values for the couplings. In

this work, we will refer to these eigenvectors as EC eigenvectors or training eigenvectors.

Once we obtain these eigenvectors, we can form them into a projection operator P and project

the Hamiltonian H at some target coupling C(cid:12) into this lower-dimensional subspace. We
call the projected Hamiltonian M = P†HP , and the norm matrix N = P†P . Then, we solve

the generalized eigenvalue problem

M(cid:126)v = λiN(cid:126)v

(3.6)

The lowest eigenvalue of this system is called the EC estimate for the ground state energy,

and the ground state eigenvector in the full space can be reconstructed from the projection

operator

|ψ0(C(cid:12))(cid:105) = P(cid:126)v0

(3.7)

From here, we briefly discuss how these techniques are implemented numerically, with a side

discussion about the computational costs. The computation associated with the eigenvector

continuation method hinges on the computation of the training eigenvectors for a few values

of the coupling c. If these can be accurately obtained, and if the overlaps between any two

EC training eigenvectors are not nearly 1, then the method will be applicable.

45

3.1 Mathematical Formalism

We start this discussion by considering a finite-dimensional linear space and a single param-

eter family of Hamiltonian matrices

H(c) = H0 + cH1

(3.8)

where H0 and H1 are Hermitian. Let |ψj(c)(cid:105) denote the eigenvectors of H(c) with corre-

sponding eigenvalues Ej(c). Since H(c) is Hermitian for real c, each Ej(c) is strictly real,

and thus there are no singularities on the real axis.

Now, we can expand Ej(c) and |ψj(c)(cid:105) as power series about c = 0.

∞(cid:88)
∞(cid:88)

n=0

Ej(c) =

|ψj(c)(cid:105) =

(n)
E
j

(0)cn/n!

|ψ

(n)
j

(0)(cid:105)cn/n!

(3.9)

(3.10)

n=0

where the superscript (n) refers to the n-th derivative. This is the approach of perturbation
theory. These series are said to converge for all |c| < |z|, where z and its complex conjugate

are the nearest singularities to c = 0 in the complex plane.

If we are interested in points |c| ≥ |z|, perturbation theory is no longer applicable. We

can instead construct an analytic extension of Eqn. (3.2) by constructing a new series about
the point c = w, where w is real and |w| < |z|.

46

Figure 3.1: Diagram of convergence region for Eqn. (3.2)

∞(cid:88)
∞(cid:88)

n=0

|ψj(c)(cid:105) =

|ψ

(n)
j

(w)(cid:105) =

|ψ

(n)
j

(w)(cid:105) (c − w)n /n!

|ψ

(n+m)
j

(0)(cid:105)wm/m!

m=0

(3.11)

(3.12)

Combining these, we can obtain

|ψj(c)(cid:105) =

∞(cid:88)

∞(cid:88)

(c − w)n wm|ψ

(n+m)
j

(0)(cid:105)/m!n!

(3.13)

n=0

m=0

that is, the eigenvector |ψj(c)(cid:105), can be expressed as a linear combination of the eigenvectors
|ψj(0)(cid:105) in the region |c| < |z| and the circular region |c − w| < |z − w| centered at w.

Analytic function theory help us understand why this technique works. Although in the

47

Figure 3.2: Diagram of Extended Convergence Region, Eqn. (3.5)

Figure 3.3: Extending the Convergence region to include the point c = v

48

Figure 3.4: Final region of convergence

original expansion for |ψj(c)(cid:105) failed to converge for |c| < |z|, using the analytic extension

Eqn. (3.5) allowed us to express a new series in terms of the original series, to arbitrary

accuracy. This process can be repeated, until the target value of the coupling c = v is within

the analytically-extended convergence region. The number of vectors required is determined

by the number of expansions needed in the analytic continuation, and the order-by-order

convergence of each series expansion.

In cases where there are singularities near the real axis, the convergence of these series

can be accelerated by including multiple eigenvectors at each value of the coupling. Let
us focus on a single eigenvector and eigenvalue of H(c), which we label as |ψ1(c)(cid:105) and E1,
respectively. Furthermore, let D be the Riemann surface associated with |ψ1(c)(cid:105). We will
assume that |ψ1(c)(cid:105) is analytic everywhere except for branch point singularities at c = z and

c = ¯z.

49

Figure 3.5: D is a finite region of the Riemann surface associated with the vector |ψ1(c)(cid:105).
We assume |ψ1(c)(cid:105) is analytic everywhere except for branch point singularities at z and ¯z.
The dashed blue line shows the path traced around z.

50

Since the matrix elements of H(c) are analytic everywhere, the characteristic polynomial

of H(c) is also analytic everywhere. We now consider what happens when traversing a closed

loop around one of these branch points, as in Fig. [3.5]. The effect this process has on the

eigenvectors and eigenvalues we will call the monodromy transformation T (z).

We consider the case in which the E1(z) is a unique eigenvalue of H(z). In this case,
T (z) transforms the eigenvector |ψ1(c)(cid:105) to a vector which is proportional to |ψ1(c)(cid:105). To
remove the proportionality constant, we pick some state |f > in the linear space such that
(cid:104)f|ψ1(c)(cid:105) (cid:54)= 0 for all points c ∈ D. This renormalized eigenvector can then be defined as

|φ1(c)(cid:105) =

1

(cid:104)f|ψ1(c)(cid:105)|ψ1(c)(cid:105)

(3.14)

Since (cid:104)f|φ1(c)(cid:105) = 1 for all c ∈ D, the transformation T (z) transforms (cid:104)f|φ1(c)(cid:105) into itself.
Therefore, |φ1(c)(cid:105) is also an invariant of the transformation. While |ψ1(c)(cid:105) has a branch

point singularity at c = z, the renormalized eigenvector now is analytic at c = z. Since
the linear space spanned by |φ1(c)(cid:105) is identical to the linear space spanned by |ψ1(c)(cid:105), it is

clear that the singularities in the eigenvector normalization have no effect on the eigenvector

continuation method.

We now consider the case in which E1(z) is not a unique eigenvalue of H(z). In this

case, the transformation T (z) can result in a permutation of eigenvalues that are degenerate

at c = z. Without loss of generality, we will assume that T (z) induces a permutation
cycle involving k eigenvalues that we label E1(c), E2(c),··· , Ek(c). With this notation, the

monodromy transformation produces the cyclic permutation

T (z) : E1(z) → E2(z) → ··· → Ek(z) → E1(z)

(3.15)

51

The case k = 1 will lead to the same conclusion as the unique eigenvalue case. We will focus
on the case k > 1. Let us label the image of |φ1(c)(cid:105) under this transformation as |φ2(c)(cid:105), an

eigenvector of H(z) with eigenvalue E2(z). We continue the labelling in this manner such

that the transformation on the eigenvectors has the form

T (z) : φ1(c)(cid:105) → φ2(c)(cid:105) → ··· → φk(c)(cid:105) → φ1(c)(cid:105)

(3.16)

The last entry in the cycle is φ1(c)(cid:105) due to the fact that (cid:104)f|φ1(c)(cid:105) is single-valued under the
transformation [T (z)]k.

We note that the linear combinations

|γn(c, z)(cid:105) =

k−1(cid:88)

j=0

ei2πnj/k|φj(k)(cid:105)

(3.17)

for n = 0,··· , k − 1 are the eigenvectors of the monodromy transformation with eigenvalues
e−i2πn/k. We can now renormalize each of these eigenvectors so that the branch point

singularity at c = z is removed for each n

|γ(cid:48)
n(c, z)(cid:105) = en log(c−z)/k|γn(c, z)(cid:105)

(3.18)

These |γ(cid:48)

n(c, z)(cid:105) vectors are now analytic at c = z. In essence, what we have done is made

linear combinations of the eigenvectors in which all of the non-integer fractional powers of
c − z have been removed. In a similar fashion, we can construct the vectors |γ(cid:48)
n(c, ¯z)(cid:105) that

are analytic at c = ¯z.

Let us now return to the discussion of performing eigenvector continuation with an EC
subspace consisting of the k eigenvectors, |ψ1(ci)(cid:105),|ψ2(ci)(cid:105),··· ,|ψk(ci)(cid:105) for each sampling

52

point ci. We note that the singularities at c = z and c = ¯z do not cause convergence
problems when when determining |ψ1(c)(cid:105). This is because |ψ1(c)(cid:105) can be written as a linear
combination of vectors |γ(cid:48)
n(c, z)(cid:105) that are analytic at c = z. Thus, to accelerate convergence

near a branch point singularity, the EC subspace should contain all the eigenvectors whose

Riemann sheets become intertwined at the branch point.

3.2 Algorithms

The strategy of eigenvector continuation is to construct a set of wavefunctions for different

values of the couplings (called the training data or training vectors,) which are used to form
a “low-dimensional” subspace that constrains the trajectory of |ψj(c)(cid:105). We start with the
smallest eigenvalue E1 and corresponding eigenvector |ψ1(c)(cid:105) of the Hamiltonian H(c). The
approach will be to sample several values of the coupling c = ci, i = 1,··· , µ, where we use

µ to denote the “order” of the EC method.

The reason for this is that in the limit of small steps between the ci’s, the derivatives in the

multi-series EC expansion can be approximated by finite differences, and the space spanned

by the eigenvectors at various couplings will be approximately identical to the derivatives of

the eigenvectors at zero coupling.

This approach gives us µ eigenvectors |ψ1(ci)(cid:105) with which to form the EC subspace. The
couplings ci are chosen to be such that |ψ1(ci)(cid:105) can be accurately computed. If it is the case
that |ψ1(ci)(cid:105) cannot be accurately computed for any ci, then eigenvector continuation will

not be possible.

We assume that c = c(cid:12) is the target value of the coupling, where we wish to determine

E1(c(cid:12)) and |ψ1(c(cid:12))(cid:105), and that direct computation of these quantities is not feasible.

53

The core algorithm of eigenvector continuation starts by forming a µ× µ norm matrix N

by taking the inner products between each of the EC wavefunctions

Ni,i(cid:48) =(cid:10)ψ1(ci(cid:48))|ψ1(ci)(cid:11)

and the matrix elements of H(c(cid:12))

Mi,i(cid:48) =(cid:10)ψ1(ci(cid:48))|H(c(cid:12))|ψ1(ci)(cid:11)

(3.19)

(3.20)

Then, we can solve the generalized eigenvalue problem, which consists of solving the equation

ˆM(cid:126)v = λ ˆN(cid:126)v

(3.21)

The smallest eigenvalue λ0 is the EC estimate for E1(c(cid:12)), and the estimate for |ψ1(c(cid:12))(cid:105) can

be constructed by the equation

|ψ1(c(cid:12))(cid:105) = P(cid:126)v1

(3.22)

where P is the projection operator, an N × µ matrix whose columns are the training eigen-

vectors.

In the Monte Carlo simulations, it is sometimes advantageous to split M into two pieces;

one containing the base interactions M0, and the other containing the term that has the

coupling being varied MEC. Then, we can compute the EC estimate for the energy for

arbitrary target coupling C(cid:12) by computing the generalized eigenvalues of

(M0 + C(cid:12)MEC) · (cid:126)v = u · N · (cid:126)v

(3.23)

54

This can be done as a post-processing step, saving time during the simulation itself.

One issue that we run into while implementing the eigenvector continuation algorithm

is in the inversion of the norm matrix required for the generalized eigenvalue problem. The

norm matrix is constructed by taking the overlaps between each of the EC eigenvectors, a

la Eqn. (3.19). If these overlaps are too close to unity, i.e.

if the EC eigenvectors are too

similar, then this norm matrix will become singular, and thus non-invertible. This provides

a strict numerical bound to this technique.

To illustrate this problem, we take as an example the Lipkin model (discussed in Section

4.4). We perform a computation to 3rd order EC, computing 3 training eigenvectors with

which we form our norm matrix. These three eigenvectors are taken to have couplings that

are integer multiples of some step size, i.e. Ci = 0∆C, 1∆C, ..., i∆C. As we vary ∆C, we

expect the overlaps between these wavefunctions to become smaller, and the norm matrix

to become less singular. In Fig [3.6], we show the 2nd eigenvalue of this norm matrix.

This second eigenvalue is a diagnostic about the invertibility of the norm matrix. It will

tend towards zero if the matrix becomes singular, and will become larger as the overlaps

between the EC eigenvectors become smaller. It is difficult to tell how small of an eigenvalue

is “too small”; in our experience, it depends on the system.

To demonstrate the effect this problem has on the extracted ground state wavefunctions

and energies, we show a plot of the EC ground state energy estimate at a target coupling

of N V /
 = 1.0, s well as the overlaps of the exact ground state wavefunction with the EC

estimate for the wavefunction. In Fig. [3.7], we show the overlaps of the two wavefunctions.
In order to plot it in a readable way, we choose to plot one minus the overlap, 1− < ψ0|ψEC >

on a log scale, so the estimate can be shown to get exponentially better as the norm matrix

becomes less singular.

55

Figure 3.6: 2nd Eigenvalue of EC norm matrix, as a function of the step size between
values of the couplings used to generate the EC eigenvectors. Here, the coupling step size is
∆C = 0.001x.

56

Figure 3.7: Overlap of EC wavefunction with exact ground state.

57

Figure 3.8: EC estimate for ground state energy vs EC coupling step size

In Fig.

[3.8], we show the estimates for the ground state energy of this system. If the

norm matrix is singular or ill-conditioned, this estimate will be impossible or very noisy.

Here, we see that the small ∆C estimates are very chaotic, and it takes a while to converge

to a steady value.

Another diagnostic that is useful in determining the stability of the algorithm is the

condition number of the norm matrix. This condition number is defined as the ratio of the

58

Figure 3.9: Here, we show the condition number of the norm matrix for 3rd order EC, in
the Lipkin model. A condition number ≈ 1 means the matrix can be inverted with no loss
of precision, whereas a large condition number indicates numerical instability during the
inversion procedure.

59

norms of the largest eigenvalue to the smallest eigenvalue

κ(N ) =

|| max λi||
|| min λi||

(3.24)

which is denoted κ. This condition number usually indicates the loss of precision that any
algorithm relying on inverting the matrix would have. κ ≈ 1 would be able to be inverted
with no loss of precision, whereas κ → ∞ would lose all information. This definition is

rather loose, but still can be helpful. In Fig.

[3.9], we show the condition number of the

same Lipkin model system, at 3rd order EC, plotted as a function of the EC coupling step

size ∆C.

Welford’s online algorithm [70] was used to extract more stable estimates for the average

and error, compared to the standard Jackknife analysis. This algorithm is an online, iterative

method in which running values for the mean ¯x and variance σ2 are kept track of, and updated

for new samples x are obtained

¯xn = ¯xn−1 +

xn − ¯xn−1

n

M2,n = M2,n−1 + (xn + ¯xn−1)(xn − ¯xn)

(3.25)

σ2
n =

M2,n
n − 1

This algorithm is well-suited to run alongside Monte Carlo simulations, since averages and

errors can be updated on the fly, rather than completely re-evaluated after each sample.

We now demonstrate how this algorithm is implemented in our Monte Carlo simulations.

This is shown diagrammatically in Figs.

[3.10] and [3.11]. First, we start by running the

code in parallel, using MPI, such that we have N processors (or cores) each running a copy of

60

the code with different initial random seeds. As we generate Monte Carlo configurations, we

group together a large number of configurations, and compute an average on each core (Step

#1). Assuming each Monte Carlo configuration is sampled from the same probability distri-

bution, these averages is guaranteed by the central limit theorem to be normally distributed.

We would then be able to compute a global average and variance/standard deviation from

the data across all cores.

In the eigenvector continuation method, we need to do some computation before we

average over the different cores, however. Here, we perform the needed eigenvalue solving,

and extract the EC estimates for the ground states. Now, we have an estimate for the energy

on each core. To begin applying the online algorithm, we will first average all of the EC

estimates on each core into a single value (Step #3). Then, Eqn. (3.25) tell us how to update

the global moving average and population variance (Step #4). After this step, we go back

to Step #1, and repeat until the average is converged and the errors are low.

In Fig.

[3.12], we show a comparison of this method with the conventional jackknife

analysis. We plot the estimates of the Coulomb energy shift in 4He, computed perturbatively

(COUL) and via direct calculation (EC), as a function of the number of Monte Carlo steps.

On the top graph, both estimates are computed via a jackknife analysis over the different

CPUs. On the bottom graph, the same perturbative result is plotted, but the EC estimate

is averaged via the online algorithm.

61

Figure 3.10: Steps #1 - #3 for the implementation of the online algorithm

62

  M(1)(core 0)N(1)(core 0)M(1)(core 1)N(1)(core 1)...M(1)(core n)N(1)(core n)500 ConfigurationsM(2)(core 0)N(2)(core 0)M(2)(core 1)N(2)(core 1)M(2)(core n)N(2)(core n)......M(500)(core 0)N(500)(core 0)M(500)(core 1)N(500)(core 1)M(500)(core n)N(500)(core n)(core 0)(core 0)(core 1)(core 1)(core n)(core n)1.Average over Configurations(core 0)2. Compute Smallest Eigenvalue.........3. Average over Cores(Standard or Jackknife)(point 1)......E0(core 1)E0(core 0)E0E0¯M¯M¯M¯N¯N¯NFigure 3.11: Step #4 for the implementation of the online algorithm

63

  For a simulation of ntotal = 10,000 configurations, collecting M and N in bins of 500, gives us 20 points for E_04. Update Average and Variance with Online Algorithm(point 1)E0¯x2=¯x1+(x2−¯x1)2¯x1=x1(point 2)E0M2,2=(x2−¯x1)∗(x2−¯x2)σ22=M2,22¯x3=¯x2+(x3−¯x2)3(point 3)E0M2,3=M2,2+(x3−¯x2)∗(x3−¯x3)σ32=M2,33¯x4=¯x3+(x4−¯x3)4(point 4)E0M2,4=M2,3+(x4−¯x3)∗(x4−¯x4)σ42=M2,44...¯xn=¯xn−1+(xn−¯xn−1)n(point n)E0M2,n=M2,n−1+(xn−¯xn−1)∗(xn−¯xn)σn2=M2,nnFigure 3.12: Plots of the perturbative Coulomb energy shift in 4He computed with a jackknife
analysis (in blue), versus the eigenvector continuation estimate for the Coulomb energy shift
(orange). On the top plot, the EC estimate was done with a jackknife analysis, and on the
bottom, it was done using the online algorithm.

64

3.3

Implementation and Computational Cost

For the most part, Monte Carlo results detailed in this work were done using a modified

version of the MCLEFT code [51], developed by B. Lu (MSU). The modifications were done

in Fortran 95, and consisted of a module file containing all the necessary subroutines.

In these codes, the bulk of the work is in generating the auxiliary field configurations,

and propagating the initial single-particle wavefunctions through Euclidean time via transfer

matrices. These single particle wavefunctions have a spin index, isospin index, and in the

modified code, an EC order index. From this we, expect the computational cost of the

eigenvector continuation code to be directly proportional to the EC order.

At the middle timestep, we insert various operators to perturbatively calculate various

quantities, like energies from higher order Chiral forces, isospin symmetry breaking terms,

Galilean invariance restoration terms, etc. These involve forming Slater determinants out

of the propagated single-particle states. Due to the size of these determinants, we only

store their logarithm. To extract the EC estimates for the energy, we use the BLAS routine

DSYGV. Since our projected matrix is a Hamiltonian matrix, we only keep the real part,

since the imaginary parts should trend to zero.

In Fig. [3.13], we show the results of profiling this code with gprof for a fixed number of

MC steps. UPDATE routines involve generation of the auxiliary fields and propagation of

the wavefunctions through time. MEASURE routines involve the computation of observables

using the wavefunctions at the center time-step. Subroutines with the suffix EC are routines

used in the eigenvector continuation method. From here, we see a rapid growth in the overall

time used by the EC method by the time that we have reached 2nd order.

This rapid scaling is in part due to the exponential growth in the number of combinations

65

Figure 3.13: gprof breakdown of the computational time spent in major subroutines in the
MCLEFT code. At 1st order, the overall time spent in the EC subroutines is about 32%,
and at 2nd order, that increases to 56%

66

Routine

1st Order EC 2nd Order EC 3rd Order EC

Update HMC
Update EC
Measure EC
Update MET
Update WAV
Measure XXX
Measure ETM
Measure IYE
Measure VQN
Measure OPE
Measure KIN
Measure PGF

Other
Total

67.79
22.74
14.93
6.96
0.31
0.26
0.22
0.18
0.15
0.03
0.12
0.06
1.12
114.63

66.37
58.29
46.09
6.76
0.27
0.25
0.21
0.17
0.13
0.02
0.12
0.02
2.62
180.08

67.52
69.96
132.53
6.85
0.33
0.25
0.21
0.17
0.18
0.03
0.12
0.02
4.90
280.08

Table 3.1: Tabulated results from gprof

you can form the eigenvectors into subspaces with. If we compute N eigenvectors, and form

a subspace out of m of them, there are

# of Choices =

(cid:18)N

(cid:19)

m

=

N !

m!(N − m)!

(3.26)

number of choices, which escalates for large N . However, if we only pick a particular set of

the N eigenvectors to form the subspace out of, the computation time spent can be reduced.

We now show how the EC algorithm scales as a function of the lattice volume. We run

the code for fixed number of Monte Carlo steps (100 in this case), a fixed number of particles

(Nf = 4), and a fixed size of the lattice in the temporal direction (Lt,in = 20,Lt,out = 20).

Here, we expect the primary term to be proportional to the lattice volume L3. We take our

fitting function to be a cubic:

T (L) = aL3 + bL2 + cL

(3.27)

67

Figure 3.14: Plot of the computational cost/time used by the EC method, as a function of
the length of the edges the lattice L. This was done for a fixed number of particles (N=4),
and a fixed number of timesteps (Ltinner = 20, Ltouter = 20).

with the y-intercept set to zero. This fit is done using a least-squares regression.

We now look at the scaling of the algorithm with respect to the number of particles Nf .

Like before, we fix the number of Monte Carlo steps (100), the spatial size of the lattice

(L = 4), and the size of the lattice in the temporal direction (Lt,in = 20,Lt,out = 20). This

time, we expect a small term proportional to N 2

f , and a larger term proportional to Nf . We

EC Order a (cubic term)
0.0729667392
0.6035009969
1.4201400218

1st
2nd
3rd

b (quadratic term)

1.4274807737
-0.7136869805
-4.4940652917

c (linear term)
-3.8672850799
-1.0474019912
4.021016857

R2

0.9972665119
0.9990421919
0.9995114091

Table 3.2: Coefficients of best polynomial fits to computational cost of the EC algorithm, as
a function of the size of the lattice L. Since there are L3 points, we expect the leading term
to be proportional to L3. We don’t expect other terms, but present them for completeness.

68

Figure 3.15: Plot of the computational cost/time used by the EC method, as a function of
the number of particles in the system. This was done for fixed lattice volume (L=4), and a
fixed number of timesteps (Ltinner = 20, Ltouter = 20).

then take our fitting function to be a quadratic:

T (Nf ) = aN 2

f + bNf

(3.28)

with the intercept set to zero as well.

69

EC Order a (quadratic term)

1st
2nd
3rd

0.0649337121
0.1325730519
0.255174513

b (linear term)
1.7797578463
4.1616152597
6.6464975649

R2

0.9999319128
0.9999292763
0.9999057979

Table 3.3: Coefficients of best polynomial fits to computational cost of the EC algorithm, as
a function of the number of particles in the system Nf . We expect a small term proportional
to L2, from the matrix elements, and a large linear term, since we keep track of single-particle
wavefunctions.

70

Chapter 4

Applications to Nuclear Systems

In this chapter, we discuss all the applications of eigenvector continuation to nuclear systems

that have been investigated so far. We start with the initial investigations into the Bose-

Hubbard model and infinite neutron matter. Then, we discuss the main project of computing

the contribution of the Coulomb interaction to the ground state energy in light nuclei. Then,

we look at a much older model of nuclear structure; the Lipkin-Meshkov-Glick Hamiltonian.

4.1 Bose-Hubbard Model

The Bose-Hubbard model in three dimensions [30] consists of a system of identical bosons

on a three-dimensional cubic lattice. The Hamiltonian has a hopping term, controlling the

nearest-neighbour hopping of each boson, a potential term proportional to U controlling the

pairwise interaction between bosons on the same site, and a chemical potential µ, and it has

the form

(cid:88)

(cid:104)n,n’(cid:105)

H = −t

a†(n’)a(n) +

(cid:88)

n

U
2

ρ(n) [ρ(n) − 1] − µ

(cid:88)

n

ρ(n)

(4.1)

where a†(n) and a(n) are the creation and annihilation operators for bosons at the lattice
site n, the first summation is over nearest-neighbour pairs (cid:104)n, n’(cid:105), and ρ(n) = a†(n)a(n) is

the density operator.

We consider a system of four particles on a 4×4×4 lattice. To illustrate why we chose this

71

Figure 4.1: Ground state energy for the Bose-Hubbard model, computed by exact diagonal-
ization and by perturbation theory.

problem, we first try using perturbation theory to find the ground state energy eigenvalue

E0 in units of the hopping parameter t. In Fig. [4.1], we show the ground state energy E0/t,

computed directly using exact diagonalization and up to 6th order in perturbation theory,

as functions of U/t. For the perturbative results, we expanded in a power series about the
point U/t = 0.0. We see immediately that there is a critical point at about U/t = −3.8,

where perturbation theory fails to converge. We note that while perturbation theory seems

to converge to the wrong value, at higher orders, it does indeed diverge at higher orders.

This divergence is due to an avoided level crossing, where the ground state eigenvalue is

merging with the 1st excited state eigenvalue.

It makes some sense that perturbation theory fails, since the physical nature of the ground

state wavefunction is changing dramatically across this critical point. It is a Bose gas for
U/t > 0, a weakly-state for −3.8 < U/t < 0, and a strongly-bound cluster for U/t < −3.8.

72

-10-8-6-4-2 0 2 4 6-5-4-3-2-1 0 1 2aE0/tU/texact energiesperturbation order 1perturbation order 2perturbation order 3perturbation order 4perturbation order 5perturbation order 6EC with 3 sampling pointssampling pointsFigure 4.2: Ground state energy of the Bose-Hubbard model, done with up to 5th order
eigenvector continuation (5 training vectors.)

Eigenvector continuation can enter the picture, once we realise that even though the

eigenvector is changing in a linear space of huge dimension, the path it traces out only

moves significantly in a small number of orthogonal directions.

To illustrate this, we consider 3rd order eigenvector continuation, in which we take three
sampling points U/t = −5.0,−1.5, 2.0. Using the ground state wavefunctions for each cou-

pling, and using the procedure from Section 2.2, we can extract the EC estimate for E0.

This result is shown also in Fig. [4.1].

At first glance, it may appear as though we got lucky. After all, we sampled values of the

coupling from each of the three ”physically-different” regions, so it seems obvious that an

expansion using linear combinations of these three wavefunctions would cover the full range

of couplings. However, we can sample couplings from entirely within one of the regions, say
from the region U/t = −2.0 to−1.6. If we now consider eigenvector continuation with up to

73

-10-8-6-4-2 0 2 4 6-5-4-3-2-1 0 1 2bE0/tU/texact energiesEC with 1 sampling pointEC with 2 sampling pointsEC with 3 sampling pointsEC with 4 sampling pointsEC with 5 sampling pointssampling pointsFigure 4.3: Computation of the lowest two energies in the Bose-Hubbard model. Here, the
lowest two eigenvectors at each training coupling were kept and used to construct the EC
subspace.

5 vectors, from the couplings U/t = −2.0,−1.9,−1.8,−1.7,−1.6, we get the results in Fig.

[4.2]. We see that even in this case, at high enough order, the EC result is able to capture
the kink in E0 at the point U/t ≈ −3.8.

In Section 3.1, we discussed how including multiple eigenvectors at each sampling point

can accelerate the convergence of the EC results. To demonstrate this, we now wish to

compute the two lowest, even-parity states of the Bose-Hubbard system. These results are

shown in Fig. [4.3]. We can now clearly see the avoided level crossing. The couplings used are
still U/t = −2.0,−1.9,−1.8,−1.7,−1.6, however, we now take the two lowest eigenvectors

at each value for the coupling, and use these to form our EC subspace. Comparing with

the results of only including one vector per sampling point, we see that the convergence has

been improved; the third order result is much closer, and the fourth/fifth order results are

74

-10-8-6-4-2 0 2 4 6-5-4-3-2-1 0 1 2E/tU/texact energiesEC with 1 sampling pointEC with 2 sampling pointsEC with 3 sampling pointsEC with 4 sampling pointsEC with 5 sampling pointssampling pointsin quite good agreement with the exact energies.

4.2 Simulations of Neutron Matter

The second system we studied was the system of pure neutron matter. The objective was

to demonstrate the use of eigenvector continuation to a full many-body Quantum Monte

Carlo simulation. We considered systems of six or fourteen neutrons at leading order (LO)

in chiral effective field theory. While there exist computationally-friendly lattice actions, as

described in [18, 19], the more demanding lattice action [33] was instead chosen to serve as

a better test for the method.

Since we are considering only neutrons, we can simplify the lattice action by reducing the

number of independent contact interactions from two to one. In the following equations, we

†
use ai (n) and a
i (n) to denote the fermion annihilation and creation operators, with spin
component i located at the lattice site n. We use the shorthand a (n) and a† (n) to denote the

column vector of nucleon components and row vector of nucleon components, respectively.
We use a spatial lattice spacing of a = 1/100 MeV−1 = 1.97 fm, and a temporal lattice
spacing of at = 1/150 MeV−1 = 1.32 fm. We also set ¯h = c = 1.

Our free, non-relativistic lattice Hamiltonian is given by

(cid:16)
(cid:17)
a†, a

=

Hfree

(cid:88)

(−1)k
2m

(cid:88)

(cid:88)

k=0,1

n

l=1,2,3

(cid:104)

(cid:16)

a† (n)

a

n + kˆl

(cid:17)

(cid:16)

+ a

(cid:17)(cid:105)

n − kˆl

(4.2)

where ˆl = ˆ1, ˆ2, ˆ3 denote the three spatial lattice unit vectors. The nucleon mass m is taken

to be 938.92 MeV.

We use the auxiliary field formalism to include the interactions between the neutrons.

75

We then express the Euclidean time-evolution operator over Lt time steps as a product of

transfer matrix operators which depend on the auxiliary field s and the pion field π0.

(cid:16)

(cid:17)

(cid:90)

U

Lt, g2
A

=

DsDπ0 exp(cid:2)−SSS (s) − Sπ0π0 (π0)(cid:3)(cid:110)

M (Lt−1) ··· M (0)(cid:111)

(4.3)

For convenience, we have explicitly listed the dependence of this operator on the square of

the axial vector coupling g2

A. The quadratic part of the auxiliary field action is

(cid:88)

n,nt

SSS (s) =

1
2

s2 (n, nt)

(4.4)

and the quadratic part of the pion field action is

Sπ0π0 (π0) =

1
2

αtm2
π

(cid:88)
(−1)k(cid:88)
(cid:88)

n,nt

π2
0 (n, nt)

+

1
2

αt

k=0,1

n,nt

l=1,2,3

(cid:88)

(cid:104)

(cid:105)
π0(n + kˆl, nt) + π0(n − kˆl, nt)

(4.5)

(4.6)

π0 (n, nt)

with the pion mass mπ is taken to be 134.98 MeV.

The normal-ordered auxiliary field transfer matrix, for a given time step nt, is given by

(cid:104)−H(nt)(a†, a, s, π0)αt

(cid:105)

M (nt) =: exp

:

(4.7)

The Hamiltonian at the time step nt is given by

H(nt)(a†, a, s, π0)αt = Hfree

(cid:16)

(cid:17)
a†, a

(cid:16)

a†, a, s

(cid:17)

(cid:16)

a†, a, π0

(cid:17)

+ S

(nt)
π0

(4.8)

(nt)
αt + S
s

76

Figure 4.4: LO chiral forces used in neutron matter simulations.
include a contact term and one-pion exchange

Included interactions

where

(nt)
s

S

(a†, a, s) =

(nt)

π0 (a†, a, π0) =

S

(cid:88)
(cid:112)−Cαt
(cid:88)

n

gAαt
2fπ

s(n, nt)a†(n)a(n)
(cid:88)

1
2

n

l=1,2,3

[π0(n + kˆl, nt) + π0(n − kˆl, nt)]a†(n)σla(n)

(4.9)

(4.10)

Here, C is the coupling of the contact interaction, and is taken to be −4.5 × 10−5 MeV−2.

This is different from the value used in [33]; however, this value was chosen to reproduce

a more realistic equation of state for neutron matter at the densities being probed in our

simulations. σl denote the Pauli spin matrices. fπ is the pion decay constant, which we take

to be 92.2 MeV. gA is the axial-vector coupling, taken as 1.29. After integrating over the

pion fields, we obtain a one-pion exchange potential that is quadratic in gA.

We performed simulations of six and fourteen neutrons in a 43 box, extracting the ground

state wavefunctions. In physical units, L = 7.9 fm, and the corresponding number densities

are 0.012 fm3 for six neutrons and 0.028 fm3 for fourteen neutrons. Our initial and final
states are taken as free Fermi gas wavefunctions, which we denote |Φ(cid:105). We use Euclidean time

77

projection to produce the ground states for given values of g2

A, which we denote |Φ, nt, g2
A(cid:105).

(cid:16)

(cid:17)|Φ(cid:105)

|Φ, nt, g2

A(cid:105) = U

nt, g2
A

(4.11)

(4.12)

(4.13)

For large nt, |Φ, nt, g2

A(cid:105) will asymptotically approach the ground state.

First, we present the results of direct calculation without eigenvector continuation. In

these calculations, we perform Monte Carlo simulations for the pion and auxiliary fields, and

calculate the ratio

r (Lt) =

(cid:104)Φ|U(cid:0)Lt, g2
(cid:1)|Φ(cid:105)
(cid:104)Φ|U(cid:0)Lt − 1, g2
(cid:1)|Φ(cid:105)

A

A

in the large Lt limit. This ratio is converted into an energy by

E(Lt) = − log{r(Lt)}/at ≈ E0 + ce−δELtat

In Fig.

[4.5], we show the direct result for the calculation of six neutrons. The red

open squares show the lattice results versus projection time, including the one-σ error bars

reflecting the stochastic errors in the simulation. From the asymptotic form in Eqn. [4.13],

the best fit is shown as a red, solid line, and the limits of the one-standard-deviation error

bars are shown as red dashed lines. In Fig. [4.6], we show the results of the direct calculation

for fourteen neutrons. The red open squares show the lattice results versus projection time,

including one-standard-deviation error bars indicating stochastic errors. The best fit is shown

as a red solid line, and the limits of the one-standard-deviation error bands are shown as

red dashed lines. Due to large sign oscillations, it is not possible to do simulations at large

projection times. This is reflected in the large uncertainties of the ground state energy

extrapolations.

78

Figure 4.5: Direct calculation of the ground state energy for six neutrons, shown as a function
of Euclidean time. The red, open squares are the lattice results, and the solid red line is
the best fit showing the exponential decay of the signal. The dashed red line is the 1σ error
band for the fit.

Figure 4.6: Direct calculation of the ground state energy for fourteen neutrons, shown as a
function of Euclidean time. The red, open squares are the lattice results, and the solid red
line is the best fit showing the exponential decay of the signal. The dashed red line is the
1σ error band for the fit.

79

 10 12 14 16 18 20 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4energy (MeV)projection time (MeV-1)direct calculationbest fiterror bands 35 40 45 50 55 60 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4energy (MeV)projection time (MeV-1)direct calculationbest fiterror bandsFor the eigenvector continuation, we compute the inner products

Ni,i(cid:48)(nt) = (cid:104)Φ, nt, c(cid:48)

i|Φ, nt, ci(cid:105)

(4.14)

for sampling values g2

A = c1, c2, c3, where c1 = 0.25, c2 = 0.60, and c3 = 0.95. We also
calculate the matrix elements of the full transfer matrix for the target value c(cid:12) = 1.292 =

1.66,

Mi,i(cid:48)(nt) = (cid:104)Φ, nt, c(cid:48)

i|U (1, c(cid:12))|Φ, nt, ci(cid:105)

(4.15)

To extract the EC energy, we solve the generalized eigenvalue problem N−1/2M N−1/2,

finding the largest eigenvalue λ(nt), and using the relation

E(nt) = − log{λ(nt)}/at ≈ E0 + ce−δE(2nt+1)at

(4.16)

In these calculations, we consider seven different choices for the EC subspace. For the
first three, we consider the 1-D subspaces formed by a single vector |Φ, nt, ci(cid:105) for i = 1, 2, 3.
The next three consist of the 2-D subspaces formed by taking two vectors |Φ, nt, ci(cid:105) and
|Φ, nt, c(cid:48)
i(cid:105) for i (cid:54)= i(cid:48). Lastly, we consider the full 3-D subspace formed from all three vectors.

The EC results for six neutrons are shown in Fig. [4.7]. The ground state energy is shown

versus projection time using sampling data g2

A = c1, c2, c3, where c1 = 0.25, c2 = 0.60, and

c3 = 0.95. We show results for one, two and three vectors. The errors are estimated using a

jackknife analysis of the Monte Carlo data. In order to decrease the sensitivity to stochastic

noise, we have used a common δE for all cases. We have also imposed a variational constraint

that adding more vectors to the EC subspace will not increase the energy. The EC results

for fourteen neutrons are shown in Fig. [4.8], with the same choices for c1, c2, c3.

80

Figure 4.7: Eigenvector continuation results for six neutrons, including up to three vectors
in the training subspace. Open squares show the EC estimates for the ground state energy
for a particular set of vectors used to construct the subspace. The lines show the best fits of
these results.

81

 13 14 15 16 17 18 19 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4energy (MeV)projection time (MeV-1)c1c2c3c2,c3c3,c1c1,c2c1,c2,c3Figure 4.8: Eigenvector continuation results for fourteen neutrons, including up to three
vectors in the training subspace. Open squares show the EC estimates for the ground state
energy for a particular set of vectors used to construct the subspace. The lines show the
best fits of these results.

82

 44 46 48 50 52 54 56 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4energy (MeV)projection time (MeV-1)c1c2c3c2,c3c3,c1c1,c2c1,c2,c3g2
A Values E0 for 6 neut. (MeV) E0 for 14 neut. (MeV)

c1
c2
c3
c2,c3
c3,c1
c1,c2
c1,c2,c3

direct calc.

ρN

13.8(1)
13.6(2)
13.6(2)
13.6(2)
13.6(2)
13.6(2)
13.6(2)

(cid:16)+3−4

(cid:17)

12
0.012 fm−3

48.9(4)
48.4(5)
48.9(6)
48.1(6)
48.9(6)
48.0(6)
48.0(6)

(cid:16)+7−15

(cid:17)

42
0.028 fm−3

Table 4.1: Summary of the neutron matter results, for 6 and 14 neutrons, compared to the
direct Monte Carlo calculations. For reference, the relevant nuclear densities are provided.

We summarize all of the results in Table [4.1]. For reference, the nuclear densities present

in each simulation are provided. These are computed by ρN = N/V , where N is the number
of neutrons in the box, and V is the volume of the box, in units of fm−3. In both cases, we

see that the estimates for the error have improved by a factor of 10 to 20.

4.3 Coulomb Interaction in 12C and 16O

The non-perturbative simulation of the Coulomb interaction in light nuclei is the main project

described in this work. For most nuclear simulations, the Coulomb interaction is handled

perturbatively. However, due to its repulsive nature, it leads to sign problems for all but the

lightest nuclei.

We want to use eigenvector continuation to circumvent this problem, by adding the

Coulomb interaction directly to the transfer matrix via auxiliary fields, and getting ground

state wavefunctions for various couplings. Then, we can build the EC estimate for the energy

shift due to the Coulomb interaction using the method described in Section 2.4.4.

83

On the lattice, the Coulomb interaction has the form

V (|n − n(cid:48)|) =

αem

max(0.5,|n − n(cid:48)|)

(4.17)

where αem = 1/137 is the fine structure constant, and |n− n(cid:48) is the distance between lattice

sites. The denominator has the max function to prevent the potential taking an infinite

value at zero distance. The value at the origin is chosen to be V (0) = 2αem. While this

approximation may seem incorrect, it doesn’t affect the p-p scattering phase shifts or any of

the low-energy nuclear observables.

To add this interaction into the transfer matrix, we define a new auxiliary field s1,

M (s, s1, t) =: exp{−Hfreeαt +

√−αemαt

s1((cid:126)n, nt)ρ((cid:126)n, nt)} :

(4.18)

where Hfree consists of the kinetic energy term, and ρ((cid:126)n, nt) is the density operator at the

(cid:126)n

lattice site ((cid:126)n, nt). The s1 auxiliary field is pulled from the distribution

(cid:88)

P = exp (−S) = exp

s1((cid:126)n, nt)V −1((cid:126)n − (cid:126)n(cid:48))s1((cid:126)n(cid:48), nt)

(4.19)

−(cid:88)

(cid:126)n,(cid:126)n(cid:48)



following the prescription in section 2.4.4.

This calculation was implemented in the MCLEFT code [51], which uses the shuttle al-

gorithm to update the auxiliary fields and propagate the wavefunctions using the transfer

matrices. Two new modules were written, eigenvectorContinuation and eigenvectorContinu-

tationAuxiliary, which were designed to compute the ground state wavefunctions for various

couplings for this s1 field. Then, the eigenvector continuation technique could be applied us-

ing the interactions in the MCLEFT code as a baseline. Also, the MCLEFT code calculates

84

Figure 4.9: Values of the Coulomb potential on the nz=0 plane. The value at the origin is
chosen to be V (0) = 2αem. This is large enough to give good results.

85

Figure 4.10: Direct calculation of Coulomb energy shifts for 12C, at C(cid:12) = 4αem. The black
dotted line shows the perturbative results, and the orange squares show the results of the
direct calculation. We see that for large projection times the sign problem begins to show
up, before we have reached a proper plateau.

the Coulomb energy perturbatively, so that was used as a benchmark for the method.

We start with 12C, computed on an L = 5 lattice, with lattice spacings of a = 1/100
MeV−1 and at = 1/1000 MeV−1. We first motivate the need for the eigenvector continu-
ation technique, by looking at the ground state energy of 12C with Coulomb included via

direct calculation. The Coulomb interaction has a noticeable sign problem in larger nuclei, so

we tune the fine structure constant to be four times its physical value, C(cid:12) = 4αem. Results

for this is shown in Fig. [4.10].

We use perturbative results as a reference value. In general, perturbative results will not

be available, or if they are, they may not be reliable. In the case of the Coulomb interaction,

86

perturbation theory works well enough to estimate the energy that it can be used as a

benchmark. These are shown with the black dotted line.

Ideally, we would compute this energy for several projection times, and show a nice

exponential decay that we could fit the true, infinite projection time ground state energy

to. However, we see the rapid onset of the Monte Carlo sign problem, which causes large

oscillations in the amplitudes at large projection time. This causes the statistical errors to

explode, making large τ extrapolations virtually impossible.

Like we did for the neutron matter system, we turn to eigenvector continuation to save

us. We can take a set of wavefunctions for small couplings, small enough that the sign

problem doesn’t bother us, and use them to extrapolate to couplings that would be otherwise

inaccessible.

We now show the results of this EC calculation for 12C, up to 3rd order. The couplings

used in the EC subspace are CEC = 0αem, 2.0αem, and 4.0αem, and the target coupling is

C(cid:12) = 4αem. This is the same target coupling used in the direct calculation, so that we can

easily compare the accuracy and stability of the results. These EC results are plotted in Fig.

[4.11].

In yellow, we show the 2nd order results, computed using the wavefunctions at coupling

CEC = 2.0αem, and 4.0αem. In green are the 3rd order results, computed with the wave-

functions at all three couplings CEC = 0αem, 2.0αem, and 4.0αem. For small projection

times, we ran into numerical instability in the norm matrix. This is due to the fact that all

the wavefunctions start in the same initial state; thus, for small projection times, the wave-

functions haven’t gotten the chance to change much. It takes until Lt = 40 (τ = 7.88 fm) for

the 2nd order results to become stable, and Lt = 50 (τ = 9.85 fm) for the 3rd order results

to become stable. We now show these two sets of data on the same plot, in Fig. [4.12], to

87

Figure 4.11: EC Calculation of Coulomb energy shifts in 12C. Again, we show the pertur-
bative results as a reference. Early on, the EC norm matrix is very singular, so data is not
present before Lt = 40. We note that the 2nd and 3rd order results are computed at the
same projection times; they are only separated on the plot for visual clarity.

88

Figure 4.12: Combined plot showing the direct results and the EC estimates of the Coulomb
energy shifts in 12C. We can immediately see the improvement gained by the EC method,
and the overall resilience to the sign problem. We note that the three points in each pack
are computed at the same projection times; they are only separated on the plot for visual
clarity.

better compare the methods. The values on this plot are summarized in Table [4.2].

In Section 3.2, we mentioned that it is often easier to break the computation of the EC

inserted operator matrix M into two pieces

M = M0 + C(cid:12)MEC

(4.20)

so that we could perform the eigenvector simulations in the post-processing stage, for arbi-

trary target coupling.

In Fig.

[4.13], we show the results of tuning the target Coulomb strength C(cid:12) = xαem

89

Timesteps

10
20
30
40
50
60
70
80
90
100

τ (fm) Perturbative Direct Calc.
-37.012(20)
1.97
-45.239(43)
3.94
5.91
-50.209(100)
-53.435(177)
7.88
-55.648(335)
9.85
-57.268(571)
11.82
13.79
-58.397(936)
-59.341(1900)
15.76
-59.814(3352)
17.73
19.70
-61.944(4708)

-37.448(51)
-45.964(52)
-50.934(59)
-54.164(63)
-56.138(77)
-57.428(72)
-58.546(84)
-59.120(99)
-59.872(105)
-60.490(99)

2nd Order EC 3rd Order EC

-
-
-

-53.982(633)
-56.002(376)
-57.553(321)
-58.716(398)
-59.788(897)
-60.601(1169)
-61.527(1316)

-
-
-
-

-56.111(95)
-57.515(45)
-58.606(130)
-59.613(474)
-60.343(680)
-61.133(686)

Table 4.2: Summary of the Coulomb results in 12C, computed via direct calculation, and by
eigenvector continuation. Dashed entries are where a singular norm matrix prevented stable
calculations.

in the range 0.0 < x < 10.0. We show the perturbative results again as a reference. Using
perturbation theory, we expect this energy correction to be linear in the coupling E0 ∝ αem.

Thus, any deviations from this would be hints of higher-order terms.

We now show a few other tests that we performed to check the validity of the simulations.

In Fig [4.14], we show the convergence of the 12C EC ground state energy estimate as a

function of the Monte Carlo step. This is a test any Monte Carlo calculation should check

first; a good result should be converged after a suitably large number of trials. We see that

here; after 7,000 steps or so, the result doesn’t vary much, and the error estimate from the

Online algorithm seems to be approaching a steady value.

We now move on to simulations of 16O. Here, the number of particles is starting to get

high enough that simulations at large projection times become even more problematic. In

Figs.

[4.15 - 4.17], we show plots like the 12C plots; the results from direct calculation of

the Coulomb energy shifts, and a result for performing the eigenvector continuation method.

We see that for large projection times, the sign problem once again kills us, though, it is

much more sever this time around. Also, the use of the eigenvector continuation has not

90

Figure 4.13: Computation of Coulomb energy shifts in 12C for arbitrary target coupling. C(cid:12)
was varied in the range 0.0αem ≤ C(cid:12) ≤ 10.0αem. Shown in the bottom figure is a zoom of
larger coupling region, showing a clear deviation from the perturbative results.

91

Figure 4.14: Plot of the EC estimate for the ground state energy of 12C, as a function of
the MC step. Shown are also the 1σ error bars at each step. These calculations were done
at a target coupling C(cid:12) = 3αem. We see that after 7000 steps, the results are more-or-less
converged, though the 3rd order results experience some numerical instability.

only stabilized the results, but has also reduced the errors rather dramatically. The values

on these plots are summarized in Table [4.3].

We now repeat the same consistency checks that we performed for the 12C simulations.

In Fig.

[4.19], we show the convergence of the 16O EC ground state energy estimate as a

function of the Monte Carlo step. It takes a bit longer to reach a stable value, but it does

appear to be converged.

92

Figure 4.15: Direct calculation of Coulomb energy shifts for 16O

Timesteps

10
20
30
40
50
60
70
80
90
100

τ (fm) Perturbative
61.406(89)
1.97
3.94
29.60(10)
4.53(10)
5.91
-14.81(12)
7.88
-29.31(12)
9.85
11.82
-39.82(13)
-47.16(13)
13.79
-52.54(15)
15.76
17.73
-56.44(15)
-59.60(15)
19.70

Direct Calc.
61.408(28)
29.579(73)
4.54(17)
-14.74(50)
-29.12(125)
-39.93(346)
-47.25(833)
-53.02(2254)
-63.55(5091)

-159.27(426912)

2nd Order EC 3rd Order EC

61.35(18)
29.51(19)
4.43(28)
-14.94(48)
-29.56(69)
-40.34(113)
-48.11(167)
-53.92(236)
-58.35(302)
-62.04(446)

-
-
-

-14.87(17)
-29.43(18)
-40.02(38)
-47.59(49)
-52.98(58)
-57.00(66)
-56.96(83)

Table 4.3: Summary of the Coulomb results in 16O, computed via direct calculation, and by
eigenvector continuation. Dashed entries are where a singular norm matrix prevented stable
calculations.

93

Figure 4.16: EC Calculation of Coulomb energy shifts in 16O. For Lt = 10 − 30, the 3rd
order results were too unstable to extract reliable data.

94

Figure 4.17: Combined plot showing the direct results and the EC estimates of the Coulomb
energy shifts in 16O.

95

Figure 4.18: Computation of Coulomb energy shifts in 16O for arbitrary target coupling.
C(cid:12) was varied in the range 0.0αem ≤ C(cid:12) ≤ 10.0αem. Shown in the bottom figure is a zoom
of larger coupling region.

96

Figure 4.19: Plot of the EC estimate for the ground state energy of 16O, as a function of the
MC step. Shown are also the 1σ error bars at each step. These calculations were done at a
target coupling C(cid:12) = 3αem.

97

4.4 EC Solutions of the Lipkin Hamiltonian

The Lipkin-Meshkov-Glick (or just Lipkin) model [50] is a model that was developed to test

many-body approximation methods and numerical techniques. It was designed to both be

exactly solvable, yet to have non-trivial behaviour.

The model consists of N fermions distributed between two levels, each assumed to be

N-fold degenerate. Each particle is assumed to be spin-1/2, and the overall wavefunction
|J, m(cid:105) is given by a total angular momentum J and projection of the angular momentum

onto the z-axis m. The interaction is written in terms of angular momentum operators

H = 
Jz +

V
2

(J 2

+ + J 2−) +

W
2

(J+J− + J−J+)

(4.21)

where Jz, J+ and J− act according to (¯h = 1)

Jz|J, m(cid:105) = m|J, m(cid:105)

J+|J, m(cid:105) =(cid:112)J(J + 1) − m(m + 1)|J, m + 1(cid:105)
J−|J, m(cid:105) =(cid:112)J(J + 1) − m(m − 1)|J, m − 1(cid:105)

(4.22)

The first term is a count of the single particle energies, while the second and third term

are what leads to interesting behaviours. The second term, or the V term, is responsible for

moving two particles from one level to the other, thus changing m by two units. The third

term, or the W term, allows swapping an upper-level particle and a lower-level particle.

This terms that allow mixing are responsible for a quantum phase transition. This feature

of the Lipkin model has been well-studied [5] [69] [10] [72] [67] [34], and has been used as a

testing ground for studying these phase transitions.

98

Figure 4.20: Basic schematic of the Lipkin model. There are N particles distributed among
two angular momentum levels, m = ±1/2, separated by some energy gap 
.

While the dimension of this problem may seem large, we note that none of the operators

in the Hamiltonian actually change J. This means the Hamiltonian matrix will be block

diagonal. There will be Jmax = N/2 blocks, one per value of J, and each block will be

dimension 2J + 1, corresponding with the values of m. In this paper, we will only consider
the (N + 1) × (N + 1) sub-matrix corresponding with J = Jmax.

The first thing that was checked was how well perturbation theory worked in this model.

For this, we used a Matlab code, developed by A. Sarkar (publication in the works). Results

for these are shown in Figs. [4.21] and [4.22].

We now demonstrate the use of the eigenvector continuation technique to this model.

Using a short Matlab code, we computed the first five eigenvectors of the Lipkin Hamiltonian,

treating V as our tunable parameter, for five different values of V; V = 0.1, 0.2, 0.3, 0.4, 0.5.

99

Figure 4.21: Results of the N=14 particle system using perturbation theory, up to 5th order.
It is clear that perturbation theory diverges quickly for couplings far from 0.

Figure 4.22: Results of the N=30 particle system using perturbation theory, up to 5th order.
Convergence here is worse than the 14 particle system.

100

Figure 4.23: EC Estimates of the five lowest energy levels of the N=14 Lipkin model. The
exact energies, obtained via diagonalizing H, are shown on the top left. At first order and
second order, the estimates are not great, but the third order estimate (and beyond) are
extremely good.

For now, we take W = 0.

In each order of EC, which we call µ, we formed a projection operator by taking the span

of the first five eigenvectors for µ of the couplings, resulting in a 5µ-dimensional operator.
Then, we projected the Hamiltonian H(V (cid:48)) into this subspace, resulting in a new matrix ˜H.

We then diagonalized ˜H; the lowest 5 eigenvalues of this matrix we call the µ-th order EC

estimates.

In Figures [4.23] and [4.25], we show the results of these calculations. We show for

reference the exact energies, shown by the solid black lines in the top-left plot. On this

top-left plot, we also show the points that were used in the construction of the EC subspace.

The other three plots on each figure show the 1st, 2nd, and 3rd order results of the EC

method. We see that while the 1st and 2nd order results are close for couplings near the

training region, they tend to not behave for other couplings. However, when we reach 3rd

101

Figure 4.24: Error between the EC estimates and the exact eigenvalues for N=14 Lipkin
model, shown on a log scale. The four lowest eigenvalues are shown, each calculated up to
5th order. 10−14 is about the machine level of precision, so better errors could be obtained
by increasing precision.

Figure 4.25: EC Estimates of the five lowest energy levels of the N=30 Lipkin model. The
exact energies, obtained via diagonalizing H, are shown on the top left. The results converge
much more slowly, and the improvements at higher orders are noticeable.

102

Figure 4.26: Error between the EC estimates and the exact eigenvalues for N=30 Lipkin
model, shown on a log scale. The four lowest eigenvalues are shown, each calculated up to
5th order. The order-by-order convergence is more pronounced with this larger number of
particles.

Figure 4.27: EC Estimates of the five lowest energy levels of the N=30 Lipkin model, with
W = 0.4. The exact energies, obtained via diagonalizing H, are shown on the top left. The
1st, 3rd, and 5th order results are shown, because convergence is overall slower than the
results from W = 0

103

Figure 4.28: Error between the EC estimates and the exact eigenvalues for N=30 Lipkin
model, with W = 0.4, shown on a log scale. The four lowest eigenvalues are shown, each
calculated up to 5th order. While the residuals are larger, the order-by-order convergence is
still present.

order, we see almost perfect agreement, down to the level of precision in the calculation. We

computed the 4th and 5th order estimates for the energies, but did not bother to plot them

because we had already reached the numerical precision threshold.

To illustrate this agreement better, we show the residuals between the EC estimates for

the energy and the exact values

(cid:0)|Ej,EC − Ej,exact|(cid:1)

R = log10

(4.23)

plotted in Figures [4.24] and [4.26]. These are organized by eigenvalue, to show the order-
by-order convergence in the results. 10−14 is as good as Matlab can do by default, so this

put a limit on how accurately our results could be. For the N = 14 system, the convergence

was already below the machine threshold by the third order calculations. For the N = 30

104

system, these errors were much larger, allowing us to see the improvements in each order of

EC.

Now, to really test the method, we turn the W term on, taking a value of W = 0.4.

This leads to a more interesting energy level spectrum, with level crossings and the like. We

still consider N = 30. In Fig. [4.27], we show the first five eigenvalues, at first, third, and

fifth order in EC. This is because the convergence isn’t as quick as before. Looking at the

residuals, shown in Fig. [4.28], it is clear that the fifth order results are not quite converged

enough.

Having demonstrated that the EC technique works, we will now push it to the strong

coupling or large N limit. Looking at the exact spectrum, we see a very tight pinch on the

spectrum around C = 1.0. This is signalling a branch point is very close to the real axis,

which would cause the convergence of the EC method to dramatically slow down (if it works

at all.)

We now apply the EC method, up to 15th order, to this system. In Figs. [4.30] and [4.31],

we show the odd EC estimates (1st order to 11th order) for the first five eigenvalues, with

the lowest eigenvalues subtracted out. On the top-left of each, we show the exact spectrum,

with the lowest eigenvalue also subtracted out, and the training points in blue circles.

We note that the EC method is not able to extrapolate past this kink in the spectrum.

While the order-by-order estimates for the eigenvalues before that point get better and

better, they all go bad by the time they reach the branch point.

Turning on the W term, we can move the location of the kink to around C = 3. Like

before, we consider up to 15th order EC, showing the results of the first five energy levels

for odd EC orders in Figs. [4.33] and [4.34], with the corresponding residuals in Figs. [4.35].

105

Figure 4.29: Plotted here are the first 20 energy levels for the N=1000 Lipkin Model. Due
to the overall energy scale, we have subtracted out the lowest energy level. We note that the
energy levels pinch together at around C = 1.0, suggesting a nearby branch point.

106

Figure 4.30: EC estimates for the lowest 5 eigenvalues of the N=1000 Lipkin model. We
show the energy with the lowest ground state eigenvalue subtracted out, to show the kink
present in the spectrum. The EC results diverge rapidly past the branch point at C = 1.

Figure 4.31: EC estimates for the lowest 5 eigenvalues of the N=1000 Lipkin model. We
show the energy with the lowest ground state eigenvalue subtracted out, to show the kink
present in the spectrum. The EC results diverge rapidly past the branch point at C = 1.

107

Figure 4.32: Residuals for the first 4 eigenvalues for the N=1000 Lipkin model. The first 15
EC orders are shown for each eigenvalue.

Figure 4.33: EC estimates for the lowest 5 eigenvalues of the N=1000 Lipkin model, with
non-zero W. We show the energy with the lowest ground state eigenvalue subtracted out,
to show the kink present in the spectrum. The EC results diverge rapidly past the branch
point at C = 3.

108

Figure 4.34: EC estimates for the lowest 5 eigenvalues of the N=1000 Lipkin model, with
non-zero W. We show the energy with the lowest ground state eigenvalue subtracted out,
to show the kink present in the spectrum. The EC results diverge rapidly past the branch
point at C = 3.

Figure 4.35: Residuals for the first 4 eigenvalues for the N=1000 Lipkin model, with non-zero
W. The first 15 EC orders are shown for each eigenvalue.

109

Chapter 5

Concluding Remarks

In this work, we have demonstrated a new technique for evaluating the extremal eigenvalues

and eigenvectors at points inaccessible to direct calculations. This technique relies on forming

a low-dimensional space out of the eigenvectors in some coupling region, and then projecting

the Hamiltonian at some target coupling into this space. Using analytic function theory, we

showed that this projected Hamiltonian can be used to determine the analytically continued

eigenvalues and eigenvectors at the target coupling.

As an application of this method, we considered four lattice systems. First, we considered

the Bose-Hubbard model, which consists of N bosons on a lattice, with a Hamiltonian

consisting of a hopping term, and pair-wise interaction, and a chemical potential. This

system has a phase transition which cannot be captured easily in perturbation theory, when

the system goes from a weakly bound Bose gas, to a strongly bound cluster state. Using

eigenvector continuation, we showed that this phase transition could be captured. The

second system studied was pure neutron matter, done at leading order in Chiral EFT for

a system of six and fourteen neutrons. In direct calculations, the sign problem prevented

accurate results from being obtained, but by the use of eigenvector continuation, we could

obtain accurate, large projection time estimates for the energies. Next, we looked at the

non-perturbative inclusion of the Coulomb interaction into our lattice calculations. Since

the Coulomb interaction is repulsive, it has a sign problem for large systems. We showed once

110

again that the direct calculations yielded poor results, and that the use of the eigenvector

continuation technique could dramatically help the determinations of the energy. Finally,

we looked at the Lipkin model, which consists of N fermions distributed into two N-fold

degenerate angular momentum levels. This levels are coupled, and allowed to mix, giving

interesting phase transitions. Like in the Bose-Hubbard model, we showed that perturbation

theory failed to capture this phase transition, and that the eigenvector continuation method

did a good job at reproducing the energy levels.

In the future, we have two primary avenues to continue this work. One is on the mathe-

matical side, looking at proving the convergence of this method, or providing error bounds

on the estimates of the eigenvalues. The other direction this work could continue is in the

applications to other systems; ones with more impact, like in QCD or in scattering calcu-

lations. The eigenvector continuation technique is very basic, so the number of research

areas where it could be applied are quite broad. We would like to continue the search for

applications, as we continue the work on the mathematical foundations of the method.

111

BIBLIOGRAPHY

112

BIBLIOGRAPHY

[1] Gert Aarts and Ion-Olimpiu Stamatescu. Stochastic quantization at finite chemical

potential. JHEP, 09:018, 2008.

[2] Andrei Alexandru, Paulo F. Bedaque, Henry Lamm, and Scott Lawrence. Finite-Density
Monte Carlo Calculations on Sign-Optimized Manifolds. Phys. Rev., D97(9):094510,
2018.

[3] Andrei Alexandru, Manfried Faber, Ivan Horvath, and Keh-Fei Liu. Lattice QCD at

finite density via a new canonical approach. Phys. Rev., D72:114513, 2005.

[4] Y. Alhassid, D. J. Dean, S. E. Koonin, G. Lang, and W. E. Ormand. Practical solution
to the monte carlo sign problem: Realistic calculations of 54Fe. Phys. Rev. Lett., 72:613–
616, Jan 1994.

[5] J. M. Arias and J. E. Garc´ıa-Ramos. Scaling Properties of the Lipkin Model at the
Critical Point. In Symmetries and Order: Algebraic Methods in Many Body Systems:
In honor of Francesco Iachello, on the occasion of his retirement New Haven, CT,
October 5-6, 2018, 2019.

[6] G. F. Bertsch, D. J. Dean, and W. Nazarewicz. Computing atomic nuclei. SciDAC

Review, 6:42, 2007.

[7] D. Bindel, J. Demmel, and M. J. Friedman. Continuation of invariant subspaces in large

bifurcation problems. SIAM J. Scientific Computing, 30:637–656, 01 2008.

[8] D. Bindel, M. Friedman, W. Govaerts, J. Hughes, and Yu.A. Kuznetsov. Numerical
computation of bifurcations in large equilibrium systems in matlab. Journal of Compu-
tational and Applied Mathematics, 261:232 – 248, 2014.

[9] Shailesh Chandrasekharan. Fermion Bag Approach to Fermion Sign Problems. Eur.

Phys. J., A49:90, 2013.

[10] Giampaolo Co’ and Stefano De Leo. Analytical and numerical analysis of the complete

Lipkin–Meshkov–Glick Hamiltonian. Int. J. Mod. Phys., E27(05):1850039, 2018.

[11] J. Cox, C. Gattringer, K. Holland, B. Scarlet, and U.-J. Wiese. Meron-cluster solution
of fermion and other sign problems. Nuclear Physics B - Proceedings Supplements, 83-
84:777 – 791, 2000. Proceedings of the XVIIth International Symposium on Lattice
Field Theory.

113

[12] M. Cristoforetti, F. Di Renzo, A. Mukherjee, and L. Scorzato. Quantum field theories

on the Lefschetz thimble. PoS, LATTICE2013:197, 2014.

[13] Marco Cristoforetti, Francesco Di Renzo, and Luigi Scorzato. New approach to the sign
problem in quantum field theories: High density QCD on a Lefschetz thimble. Phys.
Rev., D86:074506, 2012.

[14] Philippe de Forcrand and Slavo Kratochvila. Finite density QCD with a canonical

approach. Nucl. Phys. Proc. Suppl., 153:62–67, 2006. [,62(2006)].

[15] Philippe de Forcrand and Owe Philipsen. The Chiral critical line of N(f) = 2+1 QCD

at zero and non-zero baryon density. JHEP, 01:077, 2007.

[16] H. A. Dijkstra, F. W. Wubs, A. K. Cliffe, E. Doedel, I. F. Dragomirescu, B. Eckhardt,
A. Yu. Gelfgat, A. L. Hazel, V. Lucarini, A. G. Salinger, and et al. Numerical bifur-
cation methods and their application to fluid dynamics: Analysis beyond simulation.
Communications in Computational Physics, 15(1):1–45, 2014.

[17] Simon Duane, A.D. Kennedy, Brian J. Pendleton, and Duncan Roweth. Hybrid monte

carlo. Physics Letters B, 195(2):216 – 222, 1987.

[18] S. Elhatisari, E. Epelbaum, H. Krebs, T. A. L¨ahde, D. Lee, N. Li, B.-N. Lu, U.-G.
Meißner, and G. Rupak. Ab initio calculations of the isotopic dependence of nuclear
clustering. Phys. Rev. Lett., 119:222505, Dec 2017.

[19] S. Elhatisari, D. Lee, U.-G. Meißner, and G. Rupak. Nucleon-deuteron scattering using

the adiabatic projection method. The European Physical Journal A, 52(6):174, 2016.

[20] S. Elhatisari, D. Lee, G. Rupak, E. Epelbaum, H. Krebs, T. A. L¨ahde, T. Luu, and

U.-G. Meißner. Ab-initio alpha-alpha scattering. Nature, 528:111–114, 2015.

[21] S. Elhatisari, N. Li, A. Rokash, J. M. Alarc´on, D. Du, N. Klein, B.-N. Lu, U.-G. Meißner,
E. Epelbaum, H. Krebs, T. A. L¨ahde, D. Lee, and G. Rupak. Nuclear binding near a
quantum phase transition. Phys. Rev. Lett., 117:132501, Sep 2016.

[22] Evgeny Epelbaum, Hermann Krebs, Dean Lee, and Ulf-G. Meissner. Ground state
energy of dilute neutron matter at next-to-leading order in lattice chiral effective field
theory. Eur. Phys. J., A40:199–213, 2009.

[23] Evgeny Epelbaum, Hermann Krebs, Dean Lee, and Ulf-G. Meissner. Lattice chiral
effective field theory with three-body interactions at next-to-next-to-leading order. Eur.
Phys. J., A41:125–139, 2009.

114

[24] Evgeny Epelbaum, Hermann Krebs, Timo A. L¨ahde, Dean Lee, and Ulf-G. Meißner.
Dependence of the triple-alpha process on the fundamental constants of nature. Eur.
Phys. J., A49:82, 2013.

[25] Evgeny Epelbaum, Hermann Krebs, Timo A. L¨ahde, Dean Lee, Ulf-G. Meißner, and
Gautam Rupak. Ab Initio Calculation of the Spectrum and Structure of 16O. Phys.
Rev. Lett., 112(10):102501, 2014.

[26] M. Tanabashi et al. (Particle Data Group). Review of particle physics. Phys. Rev. D,

98, 2018.

[27] Z. Fodor and S. D. Katz. A New method to study lattice QCD at finite temperature

and chemical potential. Phys. Lett., B534:87–92, 2002.

[28] D. Frame, I. He, R.and Ipsen, Da. Lee, De. Lee, and E. Rrapaj. Eigenvector continuation

with subspace learning. Phys. Rev. Lett., 121(3):032501, 2018.

[29] F. Fucito, E. Marinari, G. Parisi, and C. Rebbi. A Proposal for Monte Carlo Simulations

of Fermionic Systems. Nucl. Phys., B180:369, 1981. [,586(1980)].

[30] H. A. Gersch and G. C. Knollman. Quantum cell model for bosons. Phys. Rev., 129:959–

967, Jan 1963.

[31] S. Giorgini, L. P. Pitaevskii, and S. Stringari. Theory of ultracold atomic fermi gases.

Rev. Mod. Phys., 80:1215–1274, Oct 2008.

[32] D. J. Griffiths. Introduction to Quantum Mechanics. Cambridge University Press, 2016.

[33] M. Hjorth-Jensen and collaborators. An advanced course in computational nuclear

physics. Lect. Notes Phys., 936, 2017.

[34] Yi Huang, Tongcang Li, and Zhang-qi Yin. Symmetry-breaking dynamics of the finite-
size Lipkin-Meshkov-Glick model near ground state. Phys. Rev. A, 97(1):012115, Jan
2018.

[35] J. Hubbard. Calculation of partition functions. Phys. Rev. Lett., 3:77–78, Jul 1959.

[36] F. Iachello. Analytic description of critical point nuclei in a spherical-axially deformed

shape phase transition. Phys. Rev. Lett., 87:052502, Jul 2001.

[37] Ilse C. F. Ipsen and Dean J. Lee. Determinant Approximations. arXiv e-prints, page

arXiv:1105.0437, May 2011.

[38] N. Ishii, S. Aoki, and T. Hatsuda. The Nuclear Force from Lattice QCD. Phys. Rev.

Lett., 99:022001, 2007.

115

[39] Nico Klein, Dean Lee, Weitao Liu, and Ulf-G. Meißner. Regularization Methods for

Nuclear Lattice Effective Field Theory. Phys. Lett., B747:511–516, 2015.

[40] Nico Klein, Dean Lee, and Ulf-G. Meißner. Lattice Improvement in Lattice Effective

Field Theory. Eur. Phys. J., A54(12):233, 2018.

[41] T. A. L¨ahde, T. Luu, D. Lee, U.-G. Meißner, E. Epelbaum, H. Krebs, and G. Rupak.
Nuclear lattice simulations using symmetry-sign extrapolation. The European Physical
Journal A, 51(7):92, Jul 2015.

[42] C. Lanczos. An iteration method for the solution of the eigenvalue problem of lin-
ear differential and integral operators. Journal of Research of the National Bureau of
Standards, 45, 10 1950.

[43] Kurt Langfeld and Biagio Lucini. Form the density-of-states method to finite density

quantum field theory. Acta Phys. Polon. Supp., 9:503–514, 2016.

[44] D. Lee. Ground state energy at unitarity. Phys. Rev. C, 78:024001, Aug 2008.

[45] D. Lee. Lattice simulations for few- and many-body systems. Progress in Particle and

Nuclear Physics, 63(1):117 – 154, 2009.

[46] D. Lee. From Nuclear Forces and Effective Field Theory to Nuclear Structure and

Reactions. Acta Phys. Polon., B50(3):253, 2019.

[47] D. Lee, N. Salwen, and M. Windoloski.

Introduction to stochastic error correction

methods. Physics Letters B, 502(1):329 – 337, 2001.

[48] N. Li, S. Elhatsari, E. Epelbaum, D. Lee, B.-N. Lu, and U.-G. Meissner. Neutron-proton
scattering with lattice chiral effective field theory at next-to-next-to-next-to-leading
order. Phys. Rev. C, 98:1–75, Oct 2018.

[49] Ning Li, Serdar Elhatisari, Evgeny Epelbaum, Dean Lee, Bing-Nan Lu, and Ulf-G.

Meißner. Galilean invariance restoration on the lattice. 2019.

[50] H. Lipkin, N. Meshkov, and A. Glick. Validity of many-body approximation methods

for a solvable model. Nuclear Physics, 1965.

[51] B.-N. Lu, N. Li, S. Elhatsari, E. Epelbaum, D. Lee, , and U.-G. Meissner. Essential

elements for nuclear binding.

[52] Bing-Nan Lu, Timo A. L¨ahde, Dean Lee, and Ulf-G. Meißner. Precise determination of

lattice phase shifts and mixing angles. Phys. Lett., B760:309–313, 2016.

116

[53] Timo A. L¨ahde, Evgeny Epelbaum, Hermann Krebs, Dean Lee, Ulf-G. Meißner, and
Gautam Rupak. Lattice effective field theory for nuclei from A = 4 to A = 28. 2013.
[PoSLATTICE2013,231(2014)].

[54] Timo A. L¨ahde, Evgeny Epelbaum, Hermann Krebs, Dean Lee, Ulf-G. Meißner, and
Gautam Rupak. The Hoyle state in nuclear lattice effective field theory. Pramana,
83(5):651–659, 2014.

[55] Timo A. L¨ahde, Evgeny Epelbaum, Hermann Krebs, Dean Lee, Ulf-G. Meißner, and
Gautam Rupak. Uncertainties of Euclidean Time Extrapolation in Lattice Effective
Field Theory. J. Phys., G42(3):034012, 2015.

[56] R. Machleidt and D. R. Entem. Chiral effective field theory and nuclear forces. Phys.

Rep., 503:1–75, 2011.

[57] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H.
Teller, and Edward Teller. Equation of state calculations by fast computing machines.
The Journal of Chemical Physics, 21(6):1087–1092, 1953.

[58] P.J.E. Peebles. Quantum Mechanics. Princeton University Press, 1992.

[59] M. E. Peskin and D. V. Schroeder. An Introduction to Quantum Field Theory. Avalon

Publishing, 2002.

[60] Press, Flannery, Teukolsky, and Vetterling. Numerical Recipes in C. Cambridge Uni-

versity Press, 1988.

[61] H. J. Rothe. Lattice Gauge Theories: An Introduction 4th Edition. World Scientific

Publishing Company, 2012.

[62] Y. Saad. Numerical methods for large eigenvalue problems. Classics in Applied Math-

ematics, 66, 2011.

[63] Y. Saad. Analysis of subspace iteration for eigenvalue problems with evolving matrices.

SIAM Journal on Matrix Analysis and Applications, 37:103–122, 01 2016.

[64] J.J. Sakurai and J. Napolitano. Modern Quantum Mechanics 2nd Edition. Addison-

Wesley, 2011.

[65] Stephan Scherer. Introduction to chiral perturbation theory. Adv. Nucl. Phys, 27, 2003.

[66] P. Sirkovic and D. Kressner. Subspace acceleration for large-scale parameter-dependent
hermitian eigenproblems. SIAM Journal on Matrix Analysis and Applications, 37, 04
2015.

117

[67] Pavel Str´ansk´y, Martin Dvoˇr´ak, and Pavel Cejnar. Exceptional points near first- and

second-order quantum phase transitions. Phys. Rev. E, 97(1):012112, Jan 2018.

[68] P.A. Tipler and R.A. Llewellyn. Modern Physics. W.H. Freeman and Company, 2002.

[69] Qian Wang and Francisco Perez-Bernal. Excited-state quantum phase transition and

the quantum speed limit. arXiv e-prints, page arXiv:1810.03517, Oct 2018.

[70] B. P. Welford. Note on a method for calculating corrected sums of squares and products.

Technometrics, 1926.

[71] K. G. Wilson. Confinement of quarks. Phys. Rev. D, 10:2445–2459, Oct 1974.

[72] Biao-Liang Ye, Bo Li, Xianqing Li-Jost, and Shao-Ming Fei. Quantum correlations in
critical XXZ system and LMG model. International Journal of Quantum Information,
16(3):1850029, Jan 2018.

118