ABSTRACT
Title of dissertation: ENGINEERING A QUANTUM MANY-BODY
HAMILTONIAN WITH TRAPPED IONS
Aaron Christopher Lee,
Doctor of Philosophy, 2016
Dissertation directed by: Professor Christopher Monroe
Joint Quantum Institute,
University of Maryland Department of Physics
and
National Institute of Standards and Technology
While fault-tolerant quantum computation might still be years away, analog
quantum simulators offer a way to leverage current quantum technologies to study
classically intractable quantum systems. Cutting edge quantum simulators such as
those utilizing ultracold atoms are beginning to study physics which surpass what
is classically tractable. As the system sizes of these quantum simulators increase,
there are also concurrent gains in the complexity and types of Hamiltonians which
can be simulated. In this work, I describe advances toward the realization of an
adaptable, tunable quantum simulator capable of surpassing classical computation.
We simulate long-ranged Ising and XY spin models which can have global arbitrary
transverse and longitudinal fields in addition to individual transverse fields using a
linear chain of up to 24 171Yb+ ions confined in a linear rf Paul trap. Each qubit
is encoded in the ground state hyperfine levels of an ion. Spin-spin interactions are
engineered by the application of spin-dependent forces from laser fields, coupling
spin to motion. Each spin can be read independently using state-dependent fluores-
cence. The results here add yet more tools to an ever growing quantum simulation
toolbox. One of many challenges has been the coherent manipulation of individual
qubits. By using a surprisingly large fourth-order Stark shifts in a clock-state qubit,
we demonstrate an ability to individually manipulate spins and apply independent
Hamiltonian terms, greatly increasing the range of quantum simulations which can
be implemented. As quantum systems grow beyond the capability of classical nu-
merics, a constant question is how to verify a quantum simulation. Here, I present
measurements which may provide useful metrics for large system sizes and demon-
strate them in a system of up to 24 ions during a classically intractable simulation.
The observed values are consistent with extremely large entangled states, as much
as ∼ 95% of the system entangled. Finally, we use many of these techniques in order
to generate a spin Hamiltonian which fails to thermalize during experimental time
scales due to a meta-stable state which is often called prethermal. The observed
prethermal state is a new form of prethermalization which arises due to long-range
interactions and open boundary conditions, even in the thermodynamic limit. This
prethermalization is observed in a system of up to 22 spins. We expect that system
sizes can be extended up to 30 spins with only minor upgrades to the current ap-
paratus. These results emphasize that as the technology improves, the techniques
and tools developed here can potentially be used to perform simulations which will
surpass the capability of even the most sophisticated classical techniques, enabling
the study of a whole new regime of quantum many-body physics.
ENGINEERING A QUANTUM MANY-BODY
HAMILTONIAN WITH TRAPPED IONS
by
Aaron Christopher Lee
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2016
Advisory Committee:
Professor Christopher Monroe, Chair/Advisor
Professor Andrew Childs
Professor Ian Spielman
Professor James Williams
Professor Mohammad Hafezi
© Copyright by
Aaron Christopher Lee
2016
Acknowledgments
To Ada, the love of my life
ii
Acknowledgments
When I began my journey down this road, I had no clue about what type of
people that I would meet along the way. As luck would have it, I have had the great
privilege to work beside scientists who were not only good at their job but good
people.
I first and foremost want to thank Chris Monroe, without whom none of this
would have been possible. Chris is not just a world class physicist, but also a superb
manager of his people. I have always enjoyed the benefit of his insight into field
while simultaneously able to contribute my own ideas. This is a difficult balance to
achieve, but Chris makes it look easy.
I also want to thank Phil Richerme who to this day remains one of the smartest
people I know. His positive viewpoint was a constant source of motivation which
coupled to his careful, insightful consideration of problems makes him a physicist I
emulate to this day.
Before he left us for sunny California, I had the great pleasure of working with
Wes Campbell who still stands as one of my foremost examples of what a physicist
should be. His careful, methodical approach to problems together with being one
of the most hilarious people I’ve met is something I continue to emulate today.
My years here at UMD would not have been the same without the two grad-
students who spent years putting up with my antics, Crystal Senko and Jacob Smith.
Crystal and Jake were a large part of the reason that I enjoyed going into work in
the morning, even when nothing was working. Their good humor, intelligence, and
iii
sheer determination made results which seemed impossible all too easy. My best
wishes to both of you in your endeavors.
I also want to thank the other post-docs and grad students who I have had the
privilege to work with directly: Jiehang Zhang, Paul Hess, Brian Neyenhuis, Emily
Edwards, Kihwan Kim, Simcha Korenblit, Rajibul Islam and Antonis Kyprianidis.
Each of you has taught me a lot and I thank you for your patience and guidance.
To all the other ion trappers who were always available to lend advice or
equipment, thank you for all of the help.
Without the support of our various funding agencies, this work could not have
been accomplished. Thanks to the ARO Atomic Physics Program, the AFOSR
MURI on Quantum Measurement and Verification, and the NSF Physics Frontier
Center at JQI who have supported this research.
Despite their assertions that they do not understand what I do, my parents,
Larry and Janice Lee, have always supported me throughout grad school, even
when it took their grandchildren 3000 miles away. They have provided help and
encouragement when it was needed and have made a great many sacrifices to get
me here today. My brother Alex Lee has always been a source of positive energy
these last few years. Thank you all for everything you have done for me.
Finally, I want to thank Ada, my wife, for being supportive of this endeavor
from day one, despite the sacrifices we had to make and trials we had to endure.
Without her love and support, I could not have made it through grad school. Her
insistence that I find something that I could love doing day in day out has made
my life better than I could have done on my own. On top of all of that, she is also
iv
a wonderful mother to our children: Aria, Alyssa, Andrew and Alec. They have all
been a constant well of strength and joy that I draw upon every day. Thank you
more than I can express.
v
Table of Contents
List of Figures viii
1 Introduction 1
1.1 The trapped ion quantum simulator . . . . . . . . . . . . . . . . . . . 4
1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Many-ion Hamiltonian Theory 9
2.1 Single Qubit Manipulation . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Mølmer Sørensen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Beyond the Ising Model . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1 Undesired Cross-term from TFIM . . . . . . . . . . . . . . . . 22
2.3.2 Global σz TFIM . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.3.3 Generating an X-Y Hamiltonian . . . . . . . . . . . . . . . . . 29
2.4 Trotterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.4.1 Phonon Trotterization . . . . . . . . . . . . . . . . . . . . . . 33
2.4.2 X-Y Trotterization . . . . . . . . . . . . . . . . . . . . . . . . 35
3 Experimental apparatus 39
3.1 Experimental Overview . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.1 Trapping Ions . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.2 Resonant Light . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.1.3 Longer Lifetimes . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.1.4 Coherent Control . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2 Spontaneous Emission . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.2.1 Resonant Light Leakage . . . . . . . . . . . . . . . . . . . . . 51
3.2.2 Off-resonant Spontaneous Emission . . . . . . . . . . . . . . . 53
3.2.2.1 The Two Level System . . . . . . . . . . . . . . . . . 53
3.2.2.2 Generalization to 171Yb+ . . . . . . . . . . . . . . . . 54
3.3 State Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.3.1 Camera Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.3.2 Detection Correction . . . . . . . . . . . . . . . . . . . . . . . 69
3.3.3 Independent Component Analysis (ICA) . . . . . . . . . . . . 74
vi
4 Fourth Order Stark Shift 78
4.1 Fourth-order Stark Shift Theory . . . . . . . . . . . . . . . . . . . . . 79
4.2 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.3 Experimental Demonstration . . . . . . . . . . . . . . . . . . . . . . . 90
4.4 Applications of the Fourth-order Stark Shift . . . . . . . . . . . . . . 95
5 Entanglement in Large N System 101
5.1 Spin Squeezing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.2 Quantum Fisher Information . . . . . . . . . . . . . . . . . . . . . . . 106
5.2.1 Directional QFI . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.2.2 Total QFI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.2.3 Measured QFI for 10 spins . . . . . . . . . . . . . . . . . . . . 111
5.2.4 Measured QFI for 24 spins . . . . . . . . . . . . . . . . . . . . 115
6 Failures of Quantum Thermalization 118
6.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
6.2 Observation of Prethermalization . . . . . . . . . . . . . . . . . . . . 123
6.3 Spin-boson mapping and the GGE . . . . . . . . . . . . . . . . . . . 128
6.4 Origin of the double well . . . . . . . . . . . . . . . . . . . . . . . . . 131
Bibliography 133
vii
List of Figures
2.1 Schematic representation of 171Yb+. . . . . . . . . . . . . . . . . . . . 10
2.2 Phase space evolution of a single ion. . . . . . . . . . . . . . . . . . . 16
2.3 Phase accrued from Jij. . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Normal mode contributions to Jij. . . . . . . . . . . . . . . . . . . . . 20
2.5 Trottered Spin Phase Phonon Evolution. . . . . . . . . . . . . . . . . 34
2.6 Trottered XY Spin Phase Phonon Evolution. . . . . . . . . . . . . . . 37
3.1 Schematic of trap electrodes. . . . . . . . . . . . . . . . . . . . . . . . 40
3.2 Picture of UHV chamber on table. . . . . . . . . . . . . . . . . . . . . 41
3.3 Cycling transition level diagram in 171Yb+. . . . . . . . . . . . . . . . 44
3.4 Extended diagram of 171Yb+ energy levels. . . . . . . . . . . . . . . . 47
3.5 Schematic representation of resonant light leakage diagnosis experi-
ments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.6 171Yb+ spontaneous emission level diagram. . . . . . . . . . . . . . . 55
3.7 Probability of a spontaneous emission during a pi pulse as a function
of wavelength. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.8 Probability of a single spontaneous emission in N ions. . . . . . . . . 60
3.9 Measured spontaneous emission rate. . . . . . . . . . . . . . . . . . . 61
3.10 Schematic representation of CCD hardware binning. . . . . . . . . . . 64
3.11 Example of camera summed row data. . . . . . . . . . . . . . . . . . 67
3.12 Example raw ion fidelities using ICCD. . . . . . . . . . . . . . . . . . 68
3.13 Fast ICA blind source separation. . . . . . . . . . . . . . . . . . . . . 75
3.14 Basis images from FastICA algorithm. . . . . . . . . . . . . . . . . . 76
4.1 Schematic representation of the electron energy levels of 171Yb+. . . . 81
4.2 Diagram of optics imaging 355nm light onto single ion. . . . . . . . . 83
4.3 Sketch of a typical raster pulse sequence. . . . . . . . . . . . . . . . . 88
4.4 Measured Stark shifts as a function of optical power and position. . . 91
4.5 Resolution of individual addressing beam. . . . . . . . . . . . . . . . 92
4.6 Pulse sequence for preparing a string of 10 ions in a staggered spin
configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.7 Fourth-order Stark shift used to generate disorder. . . . . . . . . . . . 96
viii
4.8 Schematic of individual rotation types. . . . . . . . . . . . . . . . . . 98
4.9 Measured Jij couplings for 7 and 10 ions. . . . . . . . . . . . . . . . . 99
5.1 Measurement of spin squeezing in 24 ions. . . . . . . . . . . . . . . . 105
5.2 Measurement of QFI in 10 Spins. . . . . . . . . . . . . . . . . . . . . 112
5.3 Measurement of QFI in 24 Spins. . . . . . . . . . . . . . . . . . . . . 116
6.1 Emergent double well potential . . . . . . . . . . . . . . . . . . . . . 121
6.2 Measured location of spin excitation . . . . . . . . . . . . . . . . . . . 123
6.3 Comparison of prethermal data to exact numerics. . . . . . . . . . . . 125
6.4 Scaling to large system size . . . . . . . . . . . . . . . . . . . . . . . 127
ix
Chapter 1: Introduction
Ever since Peter Shor introduced his now famous algorithm demonstrating
an exponential speed up in factoring large numbers [1], intense interest has been
shown in quantum information. Despite the tremendous amount of effort poured
into developing a technology capable of implementing Shor’s alogrithm and others
[2–4], fault-tolerant quantum computation remains just out of reach. Even with a
multitude of quantum error correction techniques [5], the propensity of quantum
systems toward “decoherence” remains the largest obstacle in the way of fault-
tolerance.
That is not to say that progress has not occurred. Huge strides have been
made in the engineering of quantum systems. Multiple physical implementations
exist which have the tools necessary for quantum information [6]: ultra-cold gases [7],
silicon quantum dots [8], superconducting circuits [9], and trapped ions [10], to name
a few. There are demonstrations of high fidelity single qubit gates [11–13] and two-
qubit gates [11, 12, 14] and small scale quantum algorithms are being implemented
[15,16] as the first step towards fault-tolerance. These significant milestones indicate
that fault-tolerant quantum computation is close to realization.
While fault-tolerant quantum computation will provide significant speedup
1
in special problems [4], it will also make efficient simulations of large quantum
systems possible [17]. These efficient simulations of quantum systems are critical
since a classical simulation of N simple two-level particles requires interference of
2N complex numbers, unless some approximation is used. To give some notion
of what this exponential scaling means, if we wanted to simulate a system of 275
quantum two-level particles, the amount of memory required just to store the state
of the whole system would be 7.8 × 1084 bits, assuming a double floating point
representation. Estimates of the number of atoms in the known, observable universe
range from 1078 to 1081, implying that even if every atom in the universe was used
as a classical bit of information, we still could not store the state of the system.
While there exist cases where approximations can simplify the problem, this is not
always possible, especially when large-scale entanglement is involved. A quantum
computer, on the other hand, would potentially only need to interfere N qubits,
exponentially fewer resources.
What has been described is typically known as a universal quantum computer,
yet simulation of a quantum system does not necessitate such a device. As initially
proposed by Feynman [18], a quantum simulation only requires a well-controlled
quantum system which can implement the evolution of the the desired system. This
can also be accomplished by a much simpler device which mimics the evolution in an
analog manner [19]. Due to the relative simplicity of an analog quantum simulator,
it is much easier to implement in the systems mentioned above and is a way to
surpass classical methods without requiring a universal quantum computer.
Thus, quantum simulators stand as a bridge between classical computation
2
and universal quantum computation, allowing for the computation of very specific
problems using the current cutting edge quantum technology. Quantum simulations
have been implemented in many different quantum technologies [20–23], leading to
a burgeoning field where each technology has specific types of problems that can be
explored.
One problem which is of great interest is to find the ground state of a given
Hamiltonian. If a quantum device can find the ground state of an arbitrarily con-
nected Hamiltonian in polynomial time, then such a device can be used to solve the
class of NP-complete problems [24]. Various techniques have been used to find the
ground states of specific Hamiltonians, such as quantum annealing [25,26] and adi-
abatic quantum computation [27–30]. The limitation of these techniques is that, in
general, exponentially small energy gaps are expected to prevent the global ground
state from being found.
Another significant application is to find the dynamics of a given quantum
many-body system. Graphene-like physics have been simulated with a level of con-
trol not possible in a macroscopic material [31]. Artificial gauge fields in neutral
atoms could enable the study of non-trivial topological states and even the simu-
lation of quantum chromodynamics [32]. Observations of an interacting quantum
many-body system after a quench have measured the propagation of quantum infor-
mation throughout the system [33–36]. A uniquely quantum failure of thermaliza-
tion due to many-body localization has been observed in optical lattices [37–39] and
also in trapped ions [40]. Of growing interest is the implementation of long-range
interactions, such as in dipolar molecules [41] or the van der Waals interactions in
3
Rydberg atoms [42–44]. Trapped ions also provide a platform for the study of a
long-ranged spin system with the unique ability to tune the range of the interac-
tions [45].
We are at a point in the field where systems are pushing beyond what is clas-
sically tractable while maintaining control and measurement of each qubit. In this
thesis, I will describe experiments performed on a trapped ion quantum simulator
that are on the very edge of surpassing classical techniques, having up to 24 ions
used in a classically intractable simulation. This is much larger than what has been
previously published of 18 spins [46] and can potentially even be improved in the
short term up to ∼ 30 spins.
1.1 The trapped ion quantum simulator
My focus in this thesis will be on a linear chain of trapped 171Yb+ used to sim-
ulate a many-body quantum spin Hamiltonian. We leverage two hyperfine ground
states of the ion which have a stable energy separation as our qubit or spin. Laser
beams are used to cool, prepare, and measure these states via resonant transitions
and further, apply spin-dependent forces to the ions which couple the spin to the
motion. By careful tuning of this spin-dependent force, we can only virtually couple
the spin to the motion, allowing the use of the collective motion of the atoms in the
trap as a quantum bus. This type of coupling creates a pairwise interaction graph
between the spins that falls off with ion separation r as 1/rα. The parameter α
is continuously tunable over a large domain, namely between long-range couplings,
4
which fall off slower than 1/r (small α), and short-range couplings, which fall off
faster than 1/r (large α). In addition, we can apply transverse fields in all three
directions, implement not only an Ising type Hamiltonian but also an XY Hamil-
tonain, and now even have independent control of each spin. Measurement of each
individual ion is performed using state dependent fluorescence collected with an
ICCD camera with high fidelity.
During my stay in the lab, this apparatus was used to study the ground state
frustration of an anti-ferromagnetic spin chain as a function of range [28], measure
the correlation growth rate for both an Ising and XY Hamiltonian while varying the
range across the critical range 1/r [34], apply spectroscopic techniques to verify the
interaction graph [46], produce a many-body localized state using programmable
site-specific disorder [40], and to observe a new form of prethermalization unique to
a long-ranged system with open boundaries [47]. What I think is most remarkable is
that the control of this system has significantly improved from even 3 years ago where
a 16 ion chain was the largest chain possible, to now where 24 spins is regularly used
in a classically intractable simulation while still exerting the same level of control
found in a system of 10. We think that with only a few minor improvements,
performing quantum simulations with 30 spins or more should be possible. This
system size regime is exciting as it is the boundary where exact diagonalization is
no longer possible and only numeric approximation techniques are applicable. Even
these approximations break down for certain classes of problems at > 40 spins. As I
write this, my colleagues are trapping ions in a cryogenic apparatus that should be
capable of manipulating > 50 spins with the potential to work with as many as 100
5
ions. Such a system will be more than capable of outperforming classical methods
and will represent one of the first steps beyond the limits of the classical computer.
1.2 Outline
The structure of the rest of the thesis is as follows. Chapter 2 is an in depth
discussion of how the spins are manipulated and entangled, highlighting especially
a careful derivation the interactions and their different regimes. I then discuss the
methods developed to expand the types of Hamiltonians applied, especially a tech-
nique used to produce a global σz field and how it can produce an XY interaction.
Finally I discuss some attempts to use trotterization to generate Hamiltonians and
ultimately why they didn’t work.
Chapter 3 is a very brief overview of the experimental apparatus, as extensive
documentation has been done in other theses [48–50]. I also discuss observations
about lifetimes in long chains and the current hypothesis of why what we call “de-
crystalization”, or the catastrophic loss of the ion crystal, occurs. I then derive the
spontaneous emission rate due to our Raman beams and measure the spontaneous
emission rate of ten ions to be consistent with the predicted rate. I go on to ex-
plain the current technique used to analyze the images from the ICCD camera and
the data processing method used to correct known measurement error. Finally, I
describe a method of analysis which is scalable to much larger system sizes.
Chapter 4 focuses on the surprisingly large fourth-order Stark shift imple-
mented for individual coherent control. We first lay out the underlying theory for
6
171Yb+ with a pulsed laser and then discuss the experimental implementation in
our system. We then present the observations of this fourth-order Stark shift and
demonstrate its use for state preparation and individual Hamiltonian terms.
Chapter 5 is a review of our work measuring entanglement in a large spin
system. We discuss the use of spin-squeezing as a demonstrable witness of genuine
multi-partite entanglement and measure it in a system of 24 spins, showing an
entanglement depth of 17% of the system. We then examine the quantum Fisher
information, which is a witness of genuine multi-partite entanglement assuming a
pure state, and measure this witness under the evolution of a quench to an XY
Hamiltonian which is a classically intractable problem. First a system of 10 spins
is measured and compared to numerics, where the observation is consistent with
the entire system fully entangled. This is then scaled up to a system of 24 ions
for which numerics are performed by a supercomputer, highlighting the classical
computational difficulty. Here again observations are consistent with most of the
system, up to 23 out of 24 particle entanglement.
Chapter 6 is an exposition a study of thermalization in closed quantum sys-
tems, particularly the failure of thermalization in these systems. Using the tools
from previous chapters, we are able to study both the complete failure of thermal-
ization due to many-body localization and also the suppression of thermalization
due to quasi-stationary states. I focus on the latter, which is often called prether-
malization. Previous observations of prethermal states focused on those described
by a generalized Gibbs ensemble, but here we describe and observe a new form
prethermalization that is caused by an emergent double-well potential arising from
7
the long-range interactions and open boundaries of the system, even in the ther-
modynamic limit. By changing the range of the interaction, we are able to change
the system from a prethermal state described by a generalized Gibbs ensemble to
this new kind of prethermal state. Finally, we show that this new prethermal state
persists even in large system sizes by performing the simulation on 22 spins, one of
the largest to date.
8
Chapter 2: Many-ion Hamiltonian Theory
In the following chapter, I describe the trapped ion system in the abstract and
demonstrate the theory required to generate the Hamiltonians used for quantum
simulations. I also describe some ideas explored as a means of extending control
and also the subtle reasons why some of them failed. While I describe some relevant
details here, chapter 3 will explain the experimental realization in greater detail.
The atomic system used throughout is 171Yb+ which is atomic number 70 in the
periodic table. This particular element and isotope are desirable since the nucleus
is spin-1/2 and is hydrogen-like when ionized. This makes the groundstate manifold
have a singlet/triplet structure which has a pair of “clock” states at zero magnetic
field. The qubit is encoded in these clock states, namely 2S1/2 |F = 0,mf = 0〉 (|0〉)
and 2S1/2 |F = 1,mf = 0〉 (|1〉) separated by the hyperfine frequency ωHF , because
they are magnetic field insensitive to first order which makes them ideal qubit states.
Additionally, while the resonant transition wavelengths are ultraviolet (UV) at 369
nm, they are not deep UV and can still propagate with minimal losses through air
and fibers. A more detailed discussion of the wavelengths used follows in chapter 3.
9
10
∆
ωHF2S1/2
2P1/2
2P3/2
ωZee
ωF
369 nm
355 nm
F=1
F=0
mF=-1
mF=1
mF=0
mF=0
Figure 2.1: Schematic representation of 171Yb+. The qubit is encoded in the
ground state hyperfine clock states |0〉 ≡2 S1/2 |F = 0,mf = 0〉 and |1〉 ≡2
S1/2 |F = 1,mf = 0〉. The resonant transition between the ground state 2S1/2 and
exctied state 2P1/2 manifolds is 369 nm. A pulsed laser with center frequency 355
nm is used for stimulated Raman transitions between the qubit states.
10
2.1 Single Qubit Manipulation
The foundation for any quantum information apparatus is the coherent qubit
manipulation. Coherent manipulation is what makes a pair of quantum states into
a qubit. Much of the physics described here has been discussed at length elsewhere,
but I will give a brief account again to provide a consistent notation later.
All of the coherent qubit operations are implemented using a 355 nm pulsed
laser described in more detail in section 3.1.4. Since the laser is pulsed, the optical
frequency comb is leveraged to coherently drive the qubits and even entangle them
[51]. Despite the fact that stimulated Raman transitions from two far-detuned,
non-copropagating laser beams are used for qubit manipulation rather than directly
driving the qubit frequency ωHF , the formalism is the same as that of a single beam
which has been enumerated in [52] with ~ = 1:
HI = Ω
2
σj+
(
exp
[
ı˙(η(ae−ı˙ωtt + a†eı˙ωtt)− µt+ φ)])+ h.c. (2.1)
where Ω is the Rabi frequency, a and a† are the raising and lowering operators
of the harmonic oscillator with frequency ωt, η = ∆kx0 with ∆k the difference of
the wavevectors and x0 =
√
1/2Mωt, µ = ωHF − ωL where ωL is the difference in
optical frequency of the two laser beams and φ is the difference in optical phase.
By operating in the Lamb-Dicke approximation, η  1, and the resolved sideband
11
limit, Ω ωt, Eq. 2.1 is approximated as the following:
HI ≈ Ω
2
σj+
(
1 + ı˙η(ae−ı˙ωtt + a†eı˙ωtt)
)
e−ı˙µt+ı˙φ + h.c. (2.2)
If µ = 0, the phonon terms can be approximated as zero by a rotating wave
approximation (RWA), leaving the Hamiltonian:
HCarr = Ω
2
(
σj+e
ı˙φ + σj−e
−ı˙φ) (2.3)
This interaction will flip the spin about an axis that is defined by the phase φ.
Setting φ = 0 is defined as a rotation around xˆ, whereas φ = pi/2 is a rotation
around yˆ on the Bloch sphere.
Now if µ = −ωt and RWA is applied to all terms of order ωt, the Hamiltonian
becomes
HRSB = ηΩ
2
(
aσj+e
ı˙(φ+pi/2) + a†σj−e
−ı˙(φ+pi/2)) . (2.4)
This is the standard Jaynes-Cummings Hamiltonian [53] or also known as a red
sideband interaction. If instead µ = ωt, the result is the anti-Jaynes-Cummings
Hamiltonian or blue sideband interaction, which is
HBSB = ηΩ
2
(
a†σj+e
ı˙(φ+pi/2) + aσj−e
−ı˙(φ+pi/2)) . (2.5)
12
2.2 Mølmer Sørensen
In trapped ion quantum information, the Mølmer-Sørensen (M-S) gate is the
foundation of all multi-qubit operations. This fundamental entanglement operation
has been derived many times before and in different limits. Here I will derive it again,
but the context is to clarify different approximations and explain their implications
to more sophisticated techniques.
To begin, I assume a standard Mølmer-Sørensen [54, 55] interaction arising
from the application of two off-resonant Raman transitions with the beatnote fre-
quencies detuned from both the blue sideband and red sideband. After some algebra,
the result is the following Hamiltonian [45]:
HMS =
N∑
m,j=1
ηimΩi cos [µt+ φm]σ
i
φs(ame
−iωmt + a†me
iωmt) (2.6)
where the Lamb-Dicke parameter ηi,m ≡ bi,m∆k/
√
2Mωm, bi,m is the normal mode
transformation matrix element for the ith ion and mth mode, ∆k the wave-vector
difference of along the trap axis, M the mass of a single 171Yb+, ωm is the mth
normal-mode frequency. Ωi is the ith ion two-photon Rabi frequency, {am, a†m} are
the phonon annihilation and creation operators, and σiφs = e
ı˙φsσi+ + e
−ı˙φsσi−, where
φs =
φr + φb + pi
2
φm =
φr − φb
2
(2.7)
when φr,b are the relative phases of the red and blue sideband beatnotes. Since the
13
Hamiltonian in Eq. 2.6 is time dependent, the evolution operator requires time-
ordering, which can be approximated by the Magnus expansion,
U(t) = T
[
e−ı˙
∫ t
0 dt1H(t1)
]
= eΩ1+Ω2+Ω3··· (2.8)
where T is the time-ordering operator. The first elements of the expansion are
Ω1 = −i
∫ t
0
dt1H(t1)
Ω2 = −1
2
∫ t
0
dt1
∫ t1
0
dt2[H(t1), H(t2)]
Ω3 = − i
6
∫ t
0
dt1
∫ t1
0
dt2
∫ t2
0
dt3
([
H(t1), [H(t2), H(t3)]
]
+
[
H(t3), [H(t2), H(t1)]
])
.
Calculating Ω1,
Ω1 = −
∑
i,m
σiφs(αima
†
m − α∗imam) (2.9)
where
αim =
ı˙ηimΩi
µ2 − ω2m
(
eiωmt (µ sin[µt+ φm] + ı˙ωm cos[µt+ φm])− (µ sin[φm] + iωm cos[φm])
)
.
Examining this result, the phonon operator part of Ω1 is exactly the form of a
generator of displacement. If αim is large, namely ηimΩi is large compared to µ−ωm,
then the spin and motion of the ion will be entangled. Since any measurement of the
spin traces over the phonons, when the spin and motion are entangled the observed
spin evolution appears incoherent. In fact, the trajectory of the ion in phase space
14
as a function of time and spin state can be calculated by rewriting Ω1 in terms of
the position and momentum operators,
Ω1 =
∑
i,m
σiφs
(
Re[αim](am − a†m)− ı˙Im[αim](am + a†m)
)
=ı˙
∑
i,m
σiφs
(
Re[αim]Xmpˆm − Im[αim] 1
Xm
xˆm
) (2.10)
where pˆm = −ı˙ 1Xm (am − a†m), xˆm = Xm(am + a†m), and Xm = 1√2Mωm . In Fig. 2.2,
two phase space trajectories are shown for a single ion where δ = µ − ωm has two
separate values, 3ηΩ and 0.3ηΩ. The diameter of the trajectory in phase space
is 1/δ. A small detuning thus has a larger excursion in phase space, representing
greater possible entanglement with the phonons.
Now consider the second term of the expansion, Ω2. The commutator, after
doing some algebra, is the following
[H(t1), H(t2)] =
∑
i,j,m
ηimηjmΩiΩjσ
i
φsσ
j
φs
cos (µt1 + φm) cos (µt2 + φm)
(
e−ı˙ωm(t1−t2) − eı˙ωm(t1−t2)) .
(2.11)
Performing the integration and collecting terms is a chore, but fairly straightforward.
15
-3 -2 -1 0 1 2 3
-3
-2
-1
0
1
2
3
Position
Mo
me
ntu
m 
=3,|
=3,|
=0.3,|
=0.3,|
-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3-0.3
-0.2
-0.1
0.0
0.1
0.2
0.3
Position
Mo
me
ntu
m 
=3,|
=3,|
=0.3,|
=0.3,|
Figure 2.2: Phase space evolution of a single ion. Plots of spin phase space trajec-
tories for two separate detunings, δ = 3ηΩ and δ = 0.3ηΩ where Xm ≡ 1 in the
rotating frame of trap frequency. The radius of the trajectory is ∝ 1
δ
. A non-zero fi-
nal value for xˆ and pˆ represents entanglement with the motion, with larger numbers
implying greater entanglement.
16
The final result is
Ω2 =
∑
i,j,m
σiφsσ
j
φs
ı˙ηimηjmΩiΩj
µ2 − ω2m
[
−ωmt+ ωm
2µ
sin(2φm)− ωm
2µ
sin(2(µt+ φm))
+ µ sin(φm)
(
cos((µ+ ωm)t+ φm)
µ+ ωm
− cos((µ− ωm)t+ φm)
µ− ωm
)
+ωm cos(φm)
(
sin((µ+ ωm)t+ φm)
µ+ ωm
+
sin((µ− ωm)t+ φm)
µ− ωm
)]
= ı˙
(∑
i,j
Jijσ
i
φsσ
j
φs
)
t
(2.12)
Finally, the higher order terms of the Magnus expansion are all identically zero
since all terms commute. This allows the time evolution operator to be expressed
as
U(t) = eΩ1+Ω2 . (2.13)
Ultimately, quantum information needs the time evolution operator in Eq. 2.13 to
be just the spin-spin operator in the final line of Eq. 2.12, requiring Ω1 ⇒ 0. This
can be accomplished in two different ways, resulting in two distinct regimes.
The first case is where δ = µ−ωm is the roughly the same or smaller than ηΩ.
In this “fast gate” regime, the requirement that αjm(t) = 0 is fulfilled by finding
times t such that all αjm are zero. This becomes difficult with more modes and
ions, unless one of the parameters is modulated to adjust the phase space evolution
of these motional modes. A very fruitful method is to change Ωi and the phase
φs [56, 57], which has now been used to perform arbitrary gates in a small chain of
5 ions [16].
In the second case, or the “slow gate” regime, δ is much greater than ηΩ,
17
therefore Ω1 ≈ 0. The expression for Jij can also be simplified by discarding all
terms bounded in time since δ is very large. The only remaining term is the first
in Eq. 2.12 which is linear in time. Thus, under the slow gate approximation, the
time evolution operator becomes,
U(t) = exp
(
−ı˙
∑
i,j,m
ηimηjmΩiΩj
µ2 − ω2m
ωmσ
i
φsσ
j
φs
t
)
U(t) = exp
(
−ı˙
∑
i,j,m
Jijσ
i
φsσ
j
φs
t
) (2.14)
From this expression, the effective Hamiltonian is
Heff =
∑
i,j,m
Jijσ
i
φsσ
j
φs
(2.15)
This is the foundational Hamiltonian for all trapped ion quantum simulations. Re-
lating this back to the fast-gate regime, the evolution of the modes in phase space
can be envisioned as being very small circles cycling very quickly. Since the spin-
spin interaction stems from a difference in phase between spin up and spin down
and even though the resulting circles in phase space are very small, there is a huge
number of them, producing a large phase difference and generating the effective
spin-spin coupling. Another picture is that because the interaction arises from vir-
tual coupling to phonons and no actual phonons are produced, the spins are only
virtually entangled with the motion while still generating a spin-spin interaction.
On a practical note, while a larger detuning is better for suppression of the
spin-motion entanglement, it also simultaneously suppresses Jij. Experience has
18
0 50 100 150 200
-50
-40
-30
-20
-10
0
Time (μs)
J ij
P
ha
se
0 2 4 6 8 10 12 14
-4-3
-2-1
0
δ=1.7ηΩ Approxδ=1.7ηΩδ=0.3ηΩ Approxδ=0.3ηΩ
Figure 2.3: Phase accrued from Jij. Plots of the accrued phase from Jij for two
separate detunings, δ = 1.7ηΩ and δ = 0.3ηΩ. Dashed lines are slow gate approx-
imation of the accumulated phase. Note that for short times and small detunings
(shown in inset) the approximation breaks down and very little phase is accrued.
shown that the optimal detuning for the slow gate is around 3ηΩ, giving reasonable
Jij with minimal residual spin-motion entanglement.
Since the way the two regimes are produced changes from essentially a gate
model to a time evolution model, it is often hard to connect the dots between
them. One easy way to envision the different regimes is to examine the accrual of
phase as a function of time for the full expression of Jij for different detunings. An
illustration is shown in Fig. 2.3 where the accrued phase due to Jij is shown for
two separate detunings with the corresponding slow gate approximation. Clearly
the approximation is good when δ is large, but does not work in the case of small δ
at short times. This becomes important when trying to perform dynamics at times
that are much smaller than 1/δ, as the slow-gate approximation that Jij has a purely
linear dependance on time is no longer valid.
One final observation has to do with the structure of Jij itself. The coupling
19
Mode 1
Mode 0
Mode 2
Mode 3
Mode 4
Mode 5
Mode 6
Mode 7
Mode 8
Mode 9
Jij ContributionsNormal ModeEigenvectors Cummulative Jij
Figure 2.4: Normal mode contributions to Jij. Plot of the normal mode eigenvectors
for an axial trap frequency of 0.5 MHz and radial trap frequency of 5 MHz. Mode
0 corresponds to the center-of-mass mode and mode 9 to the zig-zag mode for
a chain of ten spins. In order to visualize the contribution of each mode to the
coupling matrix, the Jij matrix contribution from each mode, ignoring the rest,
with a detuning of δ = 3ηΩ from the center-of-mass mode is shown. Also shown is
the resultant coupling matrix from a cumulative sum of the contributions. The final
coupling matrix has an effective range of α = 1.2.
20
matrix is fundamentally long-ranged, which is typically parametrized as Jij ∼ J0/|i−
j|α where |i− j| is the distance between spins i and j and α in theory ranges from
0 to 3, but in practice from about 0.5 to 1.8. The reason the interaction is long-
ranged is due to the use of the normal modes of motion as our quantum information
channel. The parameter α is varied by changing the strength of the coupling to each
mode, either by modifying the detuning or by adjusting the mode dispersion using
the axial trap frequency. An example of a short-range interaction is depicted in Fig.
2.4 showing how each mode contributes to the coupling matrix, giving the resulting
interaction graph.
2.3 Beyond the Ising Model
The easiest way to extend the simulation beyond that of a long-range Ising
model is to simply add a global transverse field. This can be accomplished by
applying a carrier that is 90 deg out of phase with the interaction, as seen from Eq.
2.3. Yet, adding a transverse field has other implications when incorporated into
the derivation of the effective Hamiltonian discussed above. Earlier, the third order
and higher terms of the Magnus expansion vanished since the commutators were
zero. With a transverse field, not only will the series no longer converge, but there
will also be an additional term in the second order of the expansion which grows
linearly with time. Below is a brief derivation to gain some insight into the size of
this lowest order unwanted term and also what can be done to mitigate it.
21
2.3.1 Undesired Cross-term from TFIM
Fixing the phase of Eq. 2.6 to φs = pi and φm = −pi/2 and also adding a
carrier term along yˆ results in the following time-dependent Hamiltonian, HTrF :
HMS =
N∑
m,j=1
ηi,mΩi sin [µt](−σjx)(ame−iωmt + a†meiωmt)
HC = ΩC
2
N∑
j
σjy
HTrF = HMS +HC .
(2.16)
When this Hamiltonian is substituted into the Magnus expansion, an additional
term in Ω1 arises, namely −ı˙
(
ΩC
2
∑N
i σ
i
y
)
t, which gives the desired transverse field
applied. There is also an unwanted term, HX , in Ω2 which comes from the spin
commutators. The size of HX can be estimated by starting with the observation,
[HTrF (t1),HTrF (t2)] = [HMS(t1),HMS(t2)] + [HMS(t1),HC(t2)]
+ [HC(t1),HMS(t2)] + [HC(t1),HC(t2)] .
(2.17)
The first commutator gives the same result as Eq. 2.12. The last commutator is
identically zero since HC is time independent. The remaining two commutators
result in an undesired cross-term, namely
HX = −1
2
∫ t
0
dt1
∫ t1
0
dt2 [HMS(t1),HC(t2)] + [HC(t1),HMS(t2)] (2.18)
22
In order to make the algebra easier, a single mode is considered and it is assumed
that δ  ηΩ, discarding all terms which have dynamics faster than δ. Thus HMS
becomes simply
HMS =
N∑
i=1
ı˙ηiΩi
2
σix(ae
ı˙δt + a†e−ı˙δt) (2.19)
After doing the commutation and integration and performing a RWA to discard all
non-stationary terms,
HX = ηiΩiΩC
4δ
σix
(
2
δ
(a− a†) + ı˙t(a+ a†)
)
(2.20)
Of course at times where t 1/δ, the only significant term is the second:
HX = ı˙
(
ηiΩiΩC
4δ
σix(a+ a
†)
)
t (2.21)
In order for HX to be negligible, Jij must fulfill the inequality Jij  ηiΩiΩC4δ which
implies that 2ηΩi  ΩC . In practice, if ηΩ ∼ 30 kHz, ΩC ≤ 20 kHz is sufficiently
small.
Since this is only the second order term and the series continues indefinitely,
it is very likely that the higher order terms have some contribution to the dynamics.
Experimental evidence suggests that these higher order terms might destructively
interfere with the second order terms since the effect of this unwanted cross-term
is observed to be smaller than expected. Even so, Eq. 2.21 still gives a reasonable
regime where HX can be neglected, and thus the final effective Hamiltonian can be
23
written
Heff =
∑
Jijσ
i
xσ
j
x +
ΩC
2
∑
i
σiy (2.22)
which is a long-ranged transverse field Ising model (TFIM).
2.3.2 Global σz TFIM
Up until the last few years, the effective Hamiltonians simulated have been
limited to an Ising Hamiltonian with either a transverse field along yˆ [28, 45] or an
additional longitudinal field [29]. The following will be a discussion of techniques
used to extend the types of Hamiltonians applied.
The first technique described will be a way to generate a significant global
σz in the system without the need to modulate the qubit frequency, since this is
difficult to do in a clock-state qubit [58]. The method itself is fairly simple. While
driving a M-S interaction, instead of having a symmetric detuning of ±µ, instead an
asymmetric detuning of ωHF −µ−D and ωHF +µ−D is used. As long as D is small
compared to δ (δ = µ− ωCM where ωCM is the center of mass frequency), then this
will result in an effective field Hz = D2
∑
i σˆ
i
z. Intuitively, one can understand this
by viewing the system from the frame defined by fixing the two sideband frequencies
in absolute frequency. In this frame of reference, if a field Bzσz was applied to the
system, relative to the sidebands, it would appear that the qubit frequency had
changed by 2Bz. This technique simply inverts the effect and generates a Bz field
by shifting the sideband detunings up or down in frequency by D, exactly as if a σz
field is applied.
24
A proper account of this technique requires that a step back from the starting
point of Eq. 2.1. Much like Eq. 2.1, the formalism will assume single beams for
brevity but can easily be extended to the stimulated Raman transitions used in the
experiment. Also assumed is that the only states coupled by the laser are the qubit
states and all other states can be neglected, namely the system can be reduced to
a two level system. Further, only one of the normal modes of motion is considered
and then generalized to many after.
Assuming a system of N qubits which have a splitting of ωHF with normal
mode frequency of ωm, the base atomic Hamiltonian is
H0 =
∑
i
ωHF
2
σiz + ωmaˆ
†
maˆm. (2.23)
Now if two lasers with frequencies ωR ≡ ωHF − ωm − δR and ωB ≡ ωHF + ωm + δB
with the same intensity are applied, the result is the following:
H =
∑
i
ωHF
2
σiz + ωmaˆ
†
maˆm +
Ω
2
σi+e
ı˙(k·r−ωRt) +
Ω
2
σi+e
ı˙(k·r−ωBt) + h.c. (2.24)
Applying the the Lamb-Dicke and resolved sideband approximations, typically the
rotating frame H0 =
ωHF
2
σiz + ωmaˆ
†
maˆm is used which, after an RWA, results in
H =
∑
i
ı˙ηΩ
2
(
σi+aˆe
ı˙δRt − σi−aˆ†e−ı˙δRt + σi+aˆ†e−ı˙δBt − σi−aˆeı˙δBt
)
(2.25)
where ωR and ωB have been substituted. If now δR = δB = δ, a standard M-S
25
Hamiltonian is recovered (see Eq. 2.19).
Instead of using the rotating frame mentioned above, now let the rotating
frame be defined H0 =
ωHF−D
2
σiz +ωmaˆ
†
maˆm. Using the same approximations as Eq.
2.25, the result is the following Hamiltonian:
H =
∑
i
D
2
σiz +
ı˙ηΩ
2
(
σi+aˆe
ı˙(δR−D)t − σi−aˆ†e−ı˙(δR−D)t
+σi+aˆ
†e−ı˙(δB+D)t − σi−aˆeı˙(δB+D)t
) (2.26)
Now the reason behind asymmetrically changing the detunings becomes clear, since
by setting δR = δ+D and δB = δ−D and substituting into Eq. 2.26 and assuming
that D  δ, the result becomes
H =
∑
i
D
2
σiz +
ı˙ηΩ
2
(
σi+aˆe
ı˙δt − σi−aˆ†e−ı˙δt + σi+aˆ†e−ı˙δt − σi−aˆeı˙δt
)
. (2.27)
The resulting Eq. 2.27 is nearly identical to Eq. 2.25 except now there is a global σz
field in addition. Reincorporating all of the motional modes, Eq. 2.27 is rewritten
in the same form as Eq. 2.6 but has an addtional global σz field. Just as discussed
in section 2.3.1, this transverse field will also generate cross-terms in the Magnus
expansion, adding a further restriction that D/2 ηΩ. When D fulfills the specified
limits, the effective Hamiltonian becomes:
Heff =
∑
Jijσ
i
xσ
j
x +
∑
i
D
2
σiz (2.28)
It should be noted here that this is a case where the rotating frame that
26
is chosen is not a “natural” frame. This technique doesn’t actually change the
hyperfine frequency, but moves the system into a frame where the qubit frequency
is different by exactly D, making the dynamics in this frame the same as a system
with the coupling matrix and a transverse field. Interestingly, this means that all
operations on the qubit, such as rotations around xˆ and yˆ on the Bloch sphere now
occur with the imposed new qubit frequency. In the lab, this means that while
normally the beat-note of the lasers is exactly ωHF in order to rotate the qubit, now
the beat-note has to be modified to ωHF −D in order to perform these rotations.
There are even more profound implications to this technique since any time
that is spent without the application of the asymmetric detunings is no longer
unitary evolution from the perspective of the effective qubit! In fact, from the
perspective of this imposed frame, what used to be unitary evolution is now the
same as a σˆz field with strength D/2. This is relevant when trying to incorporate
state preparation into the experiment since this generally takes a non-negligible
amount of time prior to the application of the interaction. If there is a delay for
some time τ prior to the interaction, the phase of the rotating frame Jij, namely
φs (see Eq. 2.7), has to be advanced by Dτ so that a σxσx interaction matches up
with the effective qubit xˆ axis and similarly for yˆ. This is equivalent to performing
a rotation around zˆ exactly equal and opposite the rotation that occurred during
the delay. If this is not done, then φs is no longer what was specified experimentally
and the observed dynamics will be incorrect. Since φs = (φR + φB + pi)/2, in order
to advance φs by Dτ , both φR and φB must increase by Dτ or just one of them
must increase by 2Dτ .
27
Similar arguments also apply to the analysis rotations which occur at the
end of the experiment where now the normal rotations around xˆ or yˆ with qubit
frequency ωHF −D also have to take into account the shifted phase from any delay
for time τ . In other words, the phase of the rotation φrot must be the same as the
total phase of the system φtot. In order for the final rotation to be coherent with
the qubit, time is tracked from the beginning of the experiment, namely for total
time T . Thus the phase of the rotation is φrot = ωRT + φ where ωR is the rotation
frequency and φ is an additional phase control. Now φtot takes into account the
evolution under each rotating frame, where if there is some delay for time τ , then
φtot = ωHF τ + (ωHF −D)(T − τ). Phase coherence requires that φrot = φtot. Since
there are two control parameters, there are two possible solutions. One possible
solution is to set ωR = ωHF −D and solve for φ:
φrot = φtot
(ωHF −D)T + φ = ωHF τ + (ωHF −D)(T − τ)
φ = Dτ
(2.29)
which is identical to what was discussed in regards to the Jij phase above. Clearly
when τ goes to zero, then this is just the same as performing the rotation with the
28
imposed frequency. The other solution is to set φ = 0 and let ωR = ωHF +D
′:
φrot = φtot
(ωHF −D′)T = ωHF τ + (ωHF −D)(T − τ)
D′T = D(T − τ)
D′ = D(1− τ
T
).
(2.30)
Both possibilities are equivalent to performing a rotation of Dτ around the zˆ axis.
These other methods are only possible because the difference between xˆ and yˆ is
just the phase of the laser beat-note as compared to the laser phase at the t = 0.
This makes the zˆ axis special since rotations can be replaced by manipulation of the
phase or frequency, over which experimental control is very good, versus performing
a rotation the lasers, which has less control due to small intensity fluctuations.
One final note is that the discussion of rotations is particularly relevant for
measurement of the xˆ and yˆ spin projections in the lab frame, yet there are no
rotations required to measure the zˆ projections in the lab frame. This is simply due
to the fact that the observable commutes with the rotating operator and thus 〈σz〉
is the same in both the lab frame and the rotating frame.
2.3.3 Generating an X-Y Hamiltonian
Here I describe an extension to a technique which produces an effective X-Y
Hamiltonian shown previously [34, 50]. What is novel here is the incorporation of
the transverse field in zˆ discussed above. Begin with the final effective Hamiltonian
29
described in Eq. 2.28 and now assume that D/2 Jij. By rewriting σixσjx in terms
of the raising and lowering operators,
Heff =
∑
Jij
(
σi+σ
j
+ + σ
i
+σ
j
− + σ
i
−σ
j
+ + σ
i
−σ
j
−
)
+
∑
i
D
2
σiz. (2.31)
If now the rotating frame of
∑
i
D
2
σiz is applied, then σ
i
± → σi±e±ı˙D/2t resulting in
Eq. 2.31 becoming
Heff =
∑
Jij
(
σi+σ
j
+e
ı˙Dt + σi+σ
j
− + σ
i
−σ
j
+ + σ
i
−σ
j
−e
−ı˙Dt) . (2.32)
Since D/2  Jij, a RW approximation of Eq. 2.32 discards both double spin-flip
terms and leaves just the spin exchange terms, which is an X-Y Hamiltonian:
Heff ≈
∑ Jij
2
(
σixσ
j
x + σ
i
yσ
j
y
)
. (2.33)
While this result is, up to a basis rotation, the same as has been previously shown,
there are unique advantages to using a transverse field in zˆ. The first is that the
strength of the field is determined by a change in beat-note frequency, which is
ultimately tied to the precision and stability of an RF source. Since these RF
sources are experimentally well controlled, the stability of this field is ultimately
limited by the noise on the qubit itself, which for a clock-state qubit is very small.
This stability makes the above approximation good out to very long times when
typical noise sources are become significant.
There is an additional experimentally relevant advantage, namely that the
30
rotations required to map the rotating frame back into the lab frame can be replaced
by ther techniques described in section 2.3.2. In fact, the rotating frame described
in Eq. 2.32 is exactly equal and opposite the frame used to produce the effective
σz field. This means that by keeping the laser beat-note at exactly ωHF for all
analysis rotations, the dynamics observed are described by Eq. 2.33 for all times.
Experimentally, this is very powerful since no calibration is required to map the
rotating frame back to the lab frame. Whereas if the transverse field is a laser
driven carrier, the calibration of the field strength is critical to the mapping back to
the lab frame and introduces noise at long times. Further, the laser driven carrier
requires measurements at only times of n2pi/B [34], restricting the observation of
the dynamics to discrete time periods.
The technique above provides an experimentally easy method of realizing an
X-Y Hamiltonian for all times, making simulations of this type of system simple to
implement even at long times. This has greatly expanded the types of simulations
which can be performed and provided a unique ability to fundamentally change the
simulated Hamiltonian.
2.4 Trotterization
Another technique explored for further expanding simulations was trotteri-
zation [17, 59]. An experimental demonstration of this technique has already been
realized in a trapped ion apparatus [60], but trotterization was revisited for two pur-
poses. The first was to enable extremely close detunings to motional modes while
31
still having negligible phonon population. The second was to have an exact X-Y
Hamiltonian without the technique explained in section 2.3.3, which by combining
the techniques would enable the production of a Heisenberg model.
The statement of trotterization is, that while for non-commuting operators A
and B, eA+B 6= eAeB, the following statement is true:
U(t) = lim
n→∞
eı˙Ht =
(
eı˙H1t/neı˙H2t/n · · · eı˙Hξt/n)n (2.34)
where H = H1 + H2 + · · · + Hξ and {H1, H2, · · · , Hξ} are a set of non-commuting
operators. Thus if n is chosen to be some large, but finite integer then
U(t) = (eı˙H1t/neı˙H2t/n · · · eı˙Hξt/n)n +∑
i<j
[Hi, Hj]
t2
2n
+
∞∑
k=3
E(k) (2.35)
where E(k) is bounded by ‖E(k)‖sup ≤ n ‖Ht/n‖ksup /k!. From this expression,
E(k) ∼ (t/n)k when k ≥ 3. If n > 1, then the sum over k in Eq. 2.35 will be smaller
than the second term which is only linearly depedent on n. Thus the dominant
source of error will be the second-order term,
∑
i<j [Hi, Hj]
t2
2n
.
This implies that to achieve any specified accuracy of the simulation of 1− ,
where  > 0, then it is required that n ∝ t2/. Whereas, solving for  and allwoing
n to change but fixing the total time t,
 ∼ t2/n = (∆t)t (2.36)
where ∆t = t/n is the time of each trotter step. The error  decreases linearly
32
with decreasing ∆t, but increases linearly with the total time of the simulation.
Essentially, for a fixed experimental time, the faster the trotterization, the better
the trotter error. As will be shown, other factors can limit the trotter step size or
even prohibit the use of trotterization entirely.
2.4.1 Phonon Trotterization
One of the first applications of trotterization stemmed from an observation that
the first term of the Magnus expansion, Ω1, is linear in σ
i
φs
while the interaction
term is quadratic. By alternating the sign of σiφs , it was supposed that any phonon
evolution would immediately be canceled by phonon evolution with opposite sign.
Trotterizing in this way would allow for detunings much closer to the motional
modes, enabling larger Jij and better control over which modes are contributing to
the coupling matrix. This intuitive picture turned out to be a naive understanding
of the problem. There are two issues which were not taken into account which make
this application not practical.
The first issue has already been discussed and illustrated in Fig. 2.3, namely
that for a small detuning δ from a particular motional mode, that mode’s sinusoidal
terms of Jij become large and cannot be ignored at times smaller than 1/δ. This
implies that the trotter time must be much larger than 1/δ, restricting how quickly
the interaction can be trottered. By itself, this fact did not prohibit the technique
since it could still be used to tune closer to the motional modes.
The second issue ends up being more important than the first. To compensate
33
-0.2 -0.1 0.0 0.1 0.2 0.3-0.5
-0.4
-0.3
-0.2
-0.1
0.0
0.1
0.2
Position
Mo
me
ntu
m 
=3
Trotter
=0.3
-0.1 0.0 0.1 0.2 0.3
-0.2
-0.1
0.0
0.1
0.2
Position
Mo
me
ntu
m 
=3
Trotter
=3
(a)
(b)
Figure 2.5: Trottered Spin Phase Phonon Evolution. Comparison of the continuous
evolution to a trottered spin phase with same detuning (a). While the phonon pop-
ulation stays small, the area enclosed in phase space, which is directly related to
the accrued Jij phase, is smaller in the trottered evolution than in the continuous
evolution. The total phase accrued by the continuous evolution is 5 times greater
than the trottered evolution. This was measured by comparing the number of ran-
domly distributed points in each shape. A similar comparison is shown (b) when the
trottered evolution has a detuning ten times smaller than the continuous evolution.
While the phonon population stays small due to the trottered phase, the enclosed
area is still 1.6 times smaller than the continuous evolution. This demonstrates that
it is better to just decrease ηΩ than to trotter the spin phase.
34
for the first issue, let the trotter step time ∆t > 1/δ. Then evolution through phase
space for the mode with the smallest detuning is no longer small circles, but star
like shapes. From the discussion of the M-S gate, generating a spin-spin interaction
is derived from a difference of phase accrued by spin-up and spin-down and so the
size of Jij is directly related to the amount of phase space enclosed by the loops in
phase space. In Fig. 2.5, the loop in phase space from a continuous evolution is
directly compared to the trotter spin phase evolution with the same detunings. Since
analytically computing the area of the arbitrary shapes is fairly difficult, instead a
uniform random distribution of points is created across the relevant area of phase
space and then the number of points within the shapes in question is counted, giving
a rough quantitative comparison of area. For the same detuning δ = 3ηΩ, the
continuous evolution is ∼ 5× larger than the trottered evolution. This is when the
detunings are the same, but even decreasing the trottered detuning to δ = 0.3ηΩ,
the continuous evolution area is still larger by a factor of ∼ 1.6×. This implies
that it would be better to decrease ηΩ by a factor of 10 and tune closer than to
trotterize the spin phase. Since it is more efficient to decrease ηΩ, there is no reason
to complicate matters by trotterizing.
2.4.2 X-Y Trotterization
The other proposed application centered around an issue that was mentioned
above, namely that a large transverse field produces undesired cross terms, even
when they are small. At long times, these cross-terms become more relevant and
35
are also phonon dependent. By trottering between a σxσx interaction and a σyσy
interaction, it was hypothesized that an X-Y Hamiltonian could be produced with-
out the cross-term or the need for rotating frames. This hypothesis was further
strengthened by the fact that an X-Y Hamiltonian had already been demonstrated
in a two-spin trapped ion qubit system using trotterization [60].
Again, the issues are more subtle than they first appear and ultimately can only
be understood in the context of the phase space picture for the motional modes. Here
the proposed experiment is to apply two separate M-S interactions, the first with
φs = pi and the second φs = pi/2, effectively trottering between a σxσx interaction
and a σyσy interaction. In the phase space picture, when φs is changed by pi/2,
|↑〉x → 1√2
(
|↑〉y + |↓〉y
)
and similarly for |↓〉x. Thus the motional population splits,
half the population maintaining the original sign and the other half swapping sign.
Because of this sign swap, the trotter time compared to 1/2δ becomes extremely
important. If the trotter time is not 1/2δ, then the phonon population stays small
and the trotterization works as expected. This can be seen in Fig. 2.6a, where
the trottered phonon population stays localized close to zero when the trotter time
τXY = .675/δ. The worst case scenario is shown in Fig. 2.6b, where the trotter time
is half 1/δ. In this case, there is some phonon population that grows linearly with
time. Each time the population splits, there are two paths in phase space. What is
computed here is only one path but it is the largest. This puts limits on the possible
size of the phonon population, but the reality is that the population fills the area
enlosed by this trajectory, which in this worst case scenario is a huge amount of
phase space. While the population going out to infinity decays exponentially with
36
-1 0 1 2 3
-3
-2
-1
0
1
2
Position
Mo
me
ntu
m

=0.31234567
8910
0 2 4 6 8
0
2
4
6
8
Position
Mo
me
ntu
m 
=0.3XY= 12 1

0.0 0.5 1.0 1.5 2.0 2.5 3.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
Position
Mo
me
ntu
m 
=0.3XY=0.675 1

(a) (b)
(c)
Figure 2.6: Trottered XY Spin Phase Phonon Evolution. Phonon population of a
single mode where the trotter time is 0.675/δ (a). This figure shows that the phonon
population is small as long as the trotter time is not close to 1/2δ. The phonon
evolution is compared to a continuous Hamiltonian with δ = 0.3ηΩ where the phonon
population is large. A similar comparison is made when the trotter time is exactly
1/(2δ) (b). Here the phonon population is growing linearly with time indicating
that certain trotter times cannot be used. Finally the phonon populations of many
modes for a random trotter time is shown. Since the populations get very large if the
trotter time is ∼ (2n+1)/(2δ) and the number of modes increases with ion number,
for any trotter time it is highly likely that some of the normal modes will have a
large phonon population. This makes trottering an XY Hamiltonian impossible for
a large system.
37
the number of trotter steps, there is still a large population entangled with the
phonons leading to effective spin decoherence.
At first glance, this issue can be solved by making the trotter time not equal
to or even close to (n + 0.5)/δ where n is an integer for all motional modes. This
is easily accomplished for two or three ions. However as more ions are added to the
trap, the density of “bad” times grows quickly due to there being more motional
modes with the same requirement and the problem rapidly becomes intractable.
An example is shown in Fig. 2.6c, where an arbitrary trotter time is chosen in a
system of ten ions. The figure shows large population excursions and are compared
to a close detuned constant interaction. Changing the time by a little changes
which modes have these large excursions, but they do not go away. For ten, it still
might be possible to find a time which does not have this phonon generation, but
any noise on the trap frequency would immediately change this stable point and
destroy the simulation. This becomes even harder as the number of spins increases.
For this reason, trotterization between a σxσx and a σyσy Hamiltonian from a M-S
interaction is not practical with a large chain of spins, despite the proof of principle
in two ions.
38
Chapter 3: Experimental apparatus
3.1 Experimental Overview
In this chapter, I will give a brief overview of the techniques and equipment
used to trap and coherently manipulate the ions. Much more detailed accounts of
this apparatus can be found in [48–50]. I will also present our current hypothesis
as to what limits the lifetimes of long chains of ions and the steps we have done in
order to mitigate it. I then show a simple derivation of the spontaneous emission
rate along with direct measurements in a chain of ten ions. Finally, I recount the
methods we use to analyzed the camera data along with an exposition of the data
processing used to recover fidelity lost by known measurement error.
3.1.1 Trapping Ions
171Yb+ ions are confined using a hand-assembled, three-layer, linear RF-Paul
trap. RF voltages with a frequency of ∼ 38 MHz are applied via a quarter-wave
helical resonator to the central layer, generating a pseudo-potential which has a
linear node along the center of the trap parallel to the electrode plane. This pseudo-
potential node defines the zˆ axis of the trap. The pseudo-potential itself is confining
39
DC
RF
DC
zx
40
0 
μm
250 μm
200 μm
Figure 3.1: Schematic of trap electrodes. RF voltage is applied to central layer
(highlighted in red) and provides confinement along xˆ and yˆ. Electrodes in gray are
grounded while electrodes in blue carry DC voltages which provide axial confinement
and also a quadrupole to break the degeneracy of the xˆ and yˆ secular frequencies.
Ions are imaged along xˆ direction.
in the radial directions (xˆ and yˆ). The outer two layers of the trap have DC voltages
applied which provide axial confinement (zˆ) and also a quadrupole used to break
the radial secular frequency degeneracy. The whole electrode assembly is suspended
in a ultra-high vacuum (UHV) chamber which is evacuated to 10−11 Torr or less.
Ions are loaded using an 171Yb isotope enriched oven which is heated to produce
a flux of neutral atoms across the trapping region. At the trapping region, neutral
atoms are excited by a laser at 399 nm, resonant with the 1S0 to
1P1 transition, which
is applied perpendicularly to the atomic flux in order to minimize Doppler shifts and
give isotope selectivity. A second beam at 355 nm photo-ionizes the atoms by way of
a resonantly assisted, dichroic, two-photon transition [61,62]. Deterministic loading
of single ions is achieved by pulsing the 355 nm light for a short time and then
40
Figure 3.2: Picture of UHV chamber on table. Chamber has a large window close to
trap for imaging ion fluorescence. Optical access is limited by microscope objective
to paths which are at ∼ 45 deg to large window via small windows on backside of
chamber. Trap axes are denoted, showing imaging axis along xˆ.
41
waiting for the ions to cool down, which is repeated until the requisite number of
ions is loaded, allowing large chains of exact size to be loaded quickly and efficiently.
By automating the loading process and using image processing techniques to count
the number of ions on the camera, I was able to greatly improve the rapidity and
ease of loading large deterministic numbers of ions.
While one ion can be loaded in a high RF trap (secular frequency in xˆ direc-
tion is ωx ∼ 5 MHz), we have found that in order to load larger chains we need to
lower both the RF confinement and the DC confinement in order to load efficiently.
We refer to this lowering of voltages as “recrystallizing”. When at the lower re-
crystallized voltages, we are able to quickly ionize and cool many ions at the cost
of generating thermal drifts in the helical resonator due to the change in applied
power. For this reason, we try to keep the time spent at the lower voltages to a
minimum and then wait for the same amount of time the trap was low after loading
in order to allow the system to re-equilibrate.
Currently what limits the number of ions which can be used in quantum sim-
ulations is collisions with the residual background gas in the UHV chamber. Under
a catastrophic collision at high RF, the chain destabilizes and “melts”. This means
that the ions are no longer able to be cooled by the lasers, which we believe is due
to the ion orbits becoming large enough that the high RF is driving energy into
the system. In order to cool the atoms and reform the linear ion crystal, we slowly
recrystallize attempting to recover all ions. For small and medium size chains, this
technique works well and allows for immediate recovery. Large chains (> 12) typi-
cally do not recool all ions. We believe that either the ions are hot enough that they
42
cannot be efficiently cooled or their spatial orbit lies outside the intensity profile of
the lasers. While it is possible that the ions not recovered might have escaped the
trap, this is unlikely due to the well depth (∼ 1 eV). Since we believe that it is more
likely that the ions are still trapped but hot, not recovering all ions requires that the
ions are dumped and a new chain loaded in order to prevent further catastrophic
events caused by the extremely hot ions in large orbits.
3.1.2 Resonant Light
The ions are cooled, initialized, and measured by means of a resonant cycling
transition between the 2S1/2 and
2P1/2 manifolds which is at 369.5 nm. With a ∼ 5
G magnetic field as a quantization axis, the Zeeman frequency is ∼ 7 MHz while
the spontaneous emission rate of 2P1/2 is 20 MHz, allowing all three Zeeman states
of the F = 1 2S1/2 to be excited using the same laser. When driven, this cycling
transition produces a high flux of photons which is then collected onto the imaging
device, whereas the singlet state of the ground state manifold cannot be coupled to
the excited state singlet because it is dipole forbidden. Thus an ion in the singlet
state appears dark, allowing for detection of the qubit state of the atom, where |1〉
is bright and |0〉 is dark.
The excited state manifold has a small probability (0.5%) of decaying into a
low lying D-state. Since the state 2D3/2 has a life-time of 52.7 ms, which is very
long compared to the experimental reprate, a repump laser at 935 nm is used to
drive population into the 3[3/2]1/2 state which has a short life-time of 37.7 ns. Since
43
12.6428  GHz2S1/2
2.105 GHz
2P1/2
γ/2π = 19.6 MHz
δZeeman = 1.4 MHz/G
36
9.
52
61
 n
m
(7
39
.0
52
1 
/ 2
)
171Yb+
F=1
F=1
F=0
F=0
2D3/2
F=2
F=1
0.86 GHz
3[3/2]1/2
F=0
F=1
2.2095 GHz
93
5.
18
79
 n
m
0
1
τ = 8.12 ns
γ/2π = 4.2 MHz
τ = 37.7 ns
γ/2π = 3.02 Hz
τ = 52.7 ms
(0.5%)
Figure 3.3: Cycling transition level diagram in 171Yb+. The primary cycling tran-
sition in 171Yb+ is between
∣∣2S1/2, F = 1〉 and ∣∣2P1/2, F = 0〉. There is a small
probability of decay into a low lying metastable state 2D3/2. Due to dipole selection
rules, this only decays to the F = 1 manifold of the D state from the F = 0 level
in the P state. A repump laser at 935 nm then drives the population back to the
F = 1 states in the ground state manifold via the F = 0 state in the 3[3/2]1/2
manifold without mixing |1〉 and |0〉 since a F = 0 to F = 0 transition is forbidden.
Since the [3/2] state has fairly high spontaneous emission rate, there is not a large
suppression of the fluorescence.
44
the coupling of the F = 0 state in the [3/2] manifold and |0〉 is dipole forbidden,
there is no mixing of |1〉 and |0〉 due to the repump, closing the cycling transition
loop. While the repump does not mix |1〉 and |0〉, there is a small probability of
off-resonant coupling via the F = 1 manifold of 2P1/2. This can lead to mis-diagnosis
of the atomic state, which will be discussed later in this chapter.
In order to initialize the ions, we resonantly couple the F = 1 manifolds of
2S1/2 and
2P1/2. The excited states all have a 1/3 probability of decaying into the
singlet ground state. Thus after scattering only a few photons, all ions are optically
pumped into |0〉 with > 99% fidelity. This enables deterministic initialization of all
spins with high fidelity. For this reason, every experiment begins by first Doppler
cooling and then optically pumping the ions into |0〉.
Doppler cooling is performed using a beam which is 10 MHz detuned from res-
onance. The optical power of the doppler cooling beam is optimized by performing
red sideband thermometry [50] on the ion during the daily calibration and adjusting
the power for minimal n¯. When holding an ion chain and not running experiments,
the Doppler cooling beam power is increased to make the chain more resilient to
background gas collisions. During loading and recrystallization, this power is fur-
ther increased to its maximum in order to quickly cool the ions and prevent ion
loss. Further, another color which is 50 MHz detuned from resonance is also applied
simultaneously in order to cool the first micro-motion sideband at 38 MHz. Without
this second color, loading and maintaining large chains becomes difficult since the
chain is not linear at the recrystallized voltages, but is in fact a zig-zag type con-
figuration. When in this configuration with large ion numbers, many ions no longer
45
sit on the RF null, making cooling more difficult due to micro-motion. Even when
the RF is restored to full strength, the chain is more resilient due to the cooling
of the first micro-motion sideband. Thus better lifetimes are achieved by cooling
both the zeroth and the first micro-motion sideband while the chain is not in use,
but during experiments only the optimized Doppler cooling that is 10 MHz detuned
is used in order to achieve the lowest possible temperature. Using this technique,
we are able to load and hold extremely large chains and still achieve extremely low
Doppler temperatures.
3.1.3 Longer Lifetimes
What follows is a small aside on observations regarding large chain lifetimes.
Recently, we were able to demonstrate the ability to coherently work with ¿20
spins. Up to now the greatest number we had been able to work with was 18 spins
[46]. The reason was that as we increased the ion number the rate of catastrophic
background gas collisions also went up. This led to a chain lifetime of ¡ 1 minute
for large crystals, whereas loading and recovering took well over a minute, making
it practically impossible to take data.
This improved dramatically once we optimized the frequency of the 935 nm
repump beam by maximizing the detection scattering rate of a chain of 10 ions.
This is very unintuitive since the 935 nm transition is power-broadened from 4 MHz
to over 150 MHz, making changes on the order of 20 MHz in the 935 frequency
imperceptible with one ion. Yet it appears that the ion chain lifetime is sensitive to
46
12.642812118466 + δ2z GHz
δ2z = (310.8)B2 Hz  [B in gauss]
2S1/2
2.105 GHz
2P1/2
γ/2π = 19.6 MHz
δZeeman = 1.4 MHz/G
36
9.
52
61
 n
m
(7
39
.0
52
1 
/ 2
)
171Yb+
F=1
F=1
F=0
F=0
2D3/2
F=2
F=1
0.86 GHz
3[3/2]1/2
F=0
F=1
2.2095 GHz
93
5.
18
79
 n
m
2F7/2
F=4
F=3
F=3
F=2
1[5/2]5/2
63
8.
61
51
 n
m
63
8.
61
02
 n
m
0
35
6.
13
4 
nm
1
2.438 µm
29
7.
1 
nm
τ = 8.12 ns
γ/2π = 4.2 MHz
τ = 37.7 ns
γ/2π = 3.02 Hz
τ = 52.7 ms
43
5.5
 nm
467 n
m
τ ~ 5.4 yrs
2D5/2 F=2
F=3
41
1.0
 nm
γ/2π = 22 Hz
τ = 7.2 ms
2P3/2
100 THz
1.35 µm
(9
9.
5%
)
(0.5%)
(0.2%
)
)
%8. 1(
(9
8.
2%
)
32
9 
nm
1.65 µm(1.0%)
(9
8.
8%
)
(17
%)
(83%
)
3.4 µm
(7/2,1)
5/2
39
7.
4 
nm
35
5 
nm
33 THz
Figure 3.4: Extended diagram of 171Yb+ energy levels. All relevant energy levels,
particularly those expected to remove ions from the cooling cycling transition. The
metastable 2F7/2 poses particular issues since once the ion is in that state it is dark
to all cooling light while having a particularly long lifetime. We believe that a
“catastrophic” collision occurs when either one or some number of ions are stuck in
this state.
the 935 frequency at the level of 20 MHz.
The current hypothesis is that the largest cause of a catastrophic collision is
that a background gas collision occurs while an ion is the 2D3/2 state, allowing for
a decay to the 2F7/2 state. This state has an extremely long life-time, preventing
that ion from undergoing Doppler cooling and effectively suppressing the cooling
rate of the chain. When the ion chain is hot, they are not well localized on the RF
null and can sample space outside the pseudo-potential approximation of the RF,
47
leading to RF heating as opposed to confinement. Due to the suppressed cooling
rate of the chain from the dark ion, chain cannot be cooled as quickly as they are
heated and so the ions are then heated until either they are too hot even for far
detuned light to cool or they are spending too much time outside the spatial profile
of the Doppler cooling beam. We believe that this is why we need to drop the RF in
order to recover or load more ions, in order that we lower the effective RF heating
rate so that our cooling can overcome the heating. This also explains why some ions
are not recovered, as they have been heated beyond our ability to cool them.
Under this hypothesis, we believe that our lifetimes are improved dramatically
by improving the repump rate simply because if there is some small probability of
collision while in the D state pc = τs where τ is the dwell time in the D state and
s is the background gas collision rate, then the probability of not having a collision
in the D state for a single ion is p = 1 − pc. Since this process is independent for
each ion, in a chain of N spins the probability that no ion has a collision is pN . This
effect makes the chain lifetimes drop off rapidly with even small decreases to this
repump rate in addition to the increased collision rate due to the increased size of
the chain. By optimizing cooling along with maximizing the repump rate, we believe
that we improved the resilience of the chain to collision, reducing the incidence of
catastrophic collisions.
48
3.1.4 Coherent Control
As mentioned in chapter 2, we use a pulsed laser to produce the necessary
frequencies used for coherent control. Specifically, we use a ND:YVO4 tripled pulsed
laser (Coherent - Paladin compact 355-4000) which has a center frequency at 355
nm. The pulse duration of the laser is τ ∼ 14 ps, giving a spectral bandwidth of
∼ 70 GHz, and has a repetition rate of frep ∼ 120.125 MHz. By selecting the comb-
tooth n = 105, we find that fHF − nfrep = 29.7 MHz. This frequency difference
is easily spanned using acousto-optical modulators (AOMs). The AOMs act as
both a fast switch and frequency control, transferring control of optical frequencies
and intensity to control of microwave frequency and power. Experimentally, this
maps the exquisite control we have in the microwave regime to the ion, allowing for
coherent control of each spin with the manipulation of microwaves.
In order to generate two beams, we split each pulse out of the laser in two using
a 50/50 beam splitter. Each beam arm is controlled using an AOM (Brimrose QZF-
210-40-355). This light is then shaped and imaged onto the ion chain where both
beams are applied in a non-copropagating configuration such that the difference in
their wave-vector (∆k) is parallel to the xˆ axis of the trap. When the difference
in the two AOMs is exactly 29.7 MHz, the two beams drive stimulated Raman
transitions. The effective Lamb-Dicke parameter of the Raman beams along the xˆ
axis is η ∼ 0.07, which is large enough to allow for a significant coupling to the red
and blue sidebands while keeping the second order sidebands minimal.
Since the 12.6 GHz beat-note at the ions is derived from the 105th comb-tooth,
49
it is extremely sensitive to small perturbations in the rep-rate, especially since we
need the beat-note to be stable down to ∼ 1 Hz. This requires a passive stability of
the rep-rate at a fraction of a Hz, which is technically challenging if not impossible.
Since the pulsed laser does not have a way to feed back to the cavity, we must
compensate for this noise in some other way. By beating a fast-photodiode signal
from the laser against a stable frequency source, we can feed forward to the frequency
of one of the Raman AOMs, locking the beatnote of the two beams [63] even as the
rep-rate fluctuates. With such a technique, the beatnote has been observed to be
stable to < 1 Hz.
3.2 Spontaneous Emission
Within quantum information, true dissipation is always a constant threat. In
trapped 171Yb+ ions, one of the primary causes of dissipation is undesired sponta-
neous emission during coherent operations. The two primary causes of this unwanted
spontaneous emission are resonant light leakage and off-resonant coupling. The first
is a problem which can always be overcome by better engineering. The second is
more difficult and due to fundamental physics and thus more interesting. In the
following, I will first give a brief account of methods which can be used to diag-
nose resonant light leakage. I will then present a derivation of spontaneous emission
due to the presence of far detuned 355 nm light and further compare the predicted
spontaneous emission rate to observed spontaneous emission in a chain of ten ions.
50
(a) DC OP τ Det
(b) DC OP R(π ) τ R(π ) Det
Figure 3.5: Schematic representation of resonant light leakage diagnosis experiments.
In this figure, “DC” represents doppler cooling, “OP” optical pumping, “Det” detec-
tion, and “R(pi)” a qubit rotation of pi. Doppler cooling sideband leakage experiment
(a). Experiment is as follows: doppler cool, optical pump, wait for time τ , detect.
Ideal measurement is zero probability bright, so any bright state population is due
to light leakage. Detection light leakage experiment (b). Experiment is as follows:
doppler cool, optical pump, pi rotation, wait for time τ , pi rotation, detect. As be-
fore, the ideal case is zero probability bright, thus bright state population is due to
leakage.
3.2.1 Resonant Light Leakage
The underlying cause of resonant light leakage is as simple as it sounds. It
means that there is light resonant with a transition that couples the qubit to an ex-
cited state, which then decays and emits a photon. The real difficulty with resonant
light leakage is that it only takes one photon from a resonant laser, and so we are
talking about optical powers on the level of less than a nW. In this case, diagnosis
is the key. There are two experiments which each test for a different kind of light
leak.
The first test is as follows (see Fig. 3.5a): Doppler cool, optical pump, wait,
and detect. The wait time τ should be on the order of 100 ms. This test measures
the spontaneous emission due to the Doppler cooling sidebands, indicating a light
leak from the Doppler cooling beam or any beam with those sidebands. The essence
of this measurement is that the sidebands will make |0〉 scatter into the F = 1
manifold of the ground state making the ion appear bright, otherwise it will remain
51
dark. Thus, if there is not any light hitting the ion, the fluorescence will remain
zero regardless of the length of τ . In practice this is fairly difficult to achieve and
so we typically try to limit it to less than 15% brightness in 100 ms.
The second test is similar (see Fig. 3.5b): Doppler cool, optical pump, rotate
qubit by pi, wait, rotate qubit by pi, detect. Again the wait time, τ , should be on
the order of 100 ms. This test measures the direct light leakage from the detection
transition by first rotating |0〉 to |1〉. If there is any detection light, the ion will
scatter with a high probability into one of the Zeeman states and so the second pi
pulse will not affect them and again the ion will appear bright. Otherwise, it will
appear dark. As before, the optimal measurement then would be zero fluorescence
regardless of wait time. It is important that the pi-pulse be well calibrated so as
not to introduce false signal, but the tolerated leakage should still be less than 15%
in 100 ms.
If any issues are observed, further identification needs to be performed to
narrow down which laser beam is leaking on the ion. The most effective way to
identify the culprit is to put a ND filter wheel into the different suspect paths to see
if the signal can be diminished or changed. Once the beam is identified, the possible
reasons for the unwanted light can be scattering off of optics, RF leakage onto the
AOM, and RF crosstalk to name a few.
52
3.2.2 Off-resonant Spontaneous Emission
Off-resonant spontaneous emission, unlike the resonant light leakage, is due
to fundamental atomic physics and as such inherent to the system. In particular, I
consider here the small off-resonant coupling due to the presence of the 355 nm light
used to generate qubit operations. Previous estimations of this rate in our system
over-estimated it by a factor of four [50, 64, 65]. I re-derive it here for clarity and
find an equation consistent with [66,67] and finally compare measured spontaneous
emission in ten ions to the predicted value.
3.2.2.1 The Two Level System
We start by assuming a simple two level system with a detuned laser applied.
The total scattering rate, γp, for this system is [68, p. 25]:
γp =
s0
γ
2
1 + s0 +
(
2∆
γ
)2 (3.1)
where γ is the linewidth of the transition, s0 ≡ I/Isat, Isat ≡ ~ω3γ12pic2 , I is the applied
intensity, ω is the radial resonant transition frequency, c is the speed of light, and ∆
is the detuning of the applied light from the resonant frequency. Now if we assume
that ∆  γ and that s0 
(
2∆
γ
)2
, namely that the detuning is much greater
than the linewidth and the saturation parameter is much smaller than the relative
53
detuning, we can make the approximation:
γp =
( γ
2∆
)2
s0
γ
2
(
1
1 +
(
γ
2∆
)2
+ s0
(
γ
2∆
)2
)
=
( γ
2∆
)2
s0
γ
2
(
1−
( γ
2∆
)2
(1 + s0) +O
( γ
2∆
)4)
= s0
γ
2
(( γ
2∆
)2
−
( γ
2∆
)4
(1 + s0) +O
( γ
2∆
)6)
≈ s0γ
2
( γ
2∆
)2
(3.2)
The single photon Rabi frequency, Γ ≡ γ
√
I
2Isat
, can be used to rewrite Eq. 3.2 as
the following:
γp = γ
(
Γ
2∆
)2
(3.3)
This gives the rate at which the far detuned light will couple to the excited state in
a simple two level system.
3.2.2.2 Generalization to 171Yb+
In 171Yb+, the laser can now couple to multiple excited states. For a first order
estimation, we will consider simply the total rate of spontaneous emission and not
worry about the final state. Under this assumption, we can calculate each excited
state independently and their total will be the total spontaneous emission rate, γp,
from a starting state A. Thus,
γ(p,A) =
∑
B
γB
(
ΓAB
2∆B
)2
(3.4)
54
01
νHF=12.6428 GHz
Δ=33.9 THz
Δ'=66.2 THz
ΔBra=16.9 THz
2S1/2
F=1
F=0
2P1/2 F=1F=0
2P3/2
F=2
F=1
3[3/2]3/2
354.7 nm
Figure 3.6: 171Yb+ spontaneous emission level diagram. Relevant electron states
when considering spontaneous emission from the 355 nm laser. Note that the
3[3/2]3/2 state has been included as it’s proximity has an effect on the spontaneous
emission of ≈ 10%.
55
where γB is the linewidth of the excited state B and ∆B is the laser’s detuning from
state B. ΓAB is the coupling between the two states by the laser:
ΓAB = 〈B| ~E · ~µ |A〉 (3.5)
where ~E = E ˆ is electric field of the laser with polarization ˆ and ~µ is the dipole
moment. Now by applying Fermi’s Golden Rule and doing some math, the following
can be derived, but this has been done by Mizrahi [65, p. 77]. I will just summarize
the results, namely
ΓAB =
(
CAB,
√
(2J ′ + 1)
)
ΓB (3.6)
where CAB, is the Clebsch-Gordan coefficients between the two states, A and B
with a laser polarization of ˆ, J ′ is the total electron angular momentum of the
excited state B, and ΓB = γB
√
I
2IsatB
is the single photon Rabi frequency for the
transition. Let us define
ˆ = −σˆ− + pipˆi + +σˆ+ (3.7)
and further require that 2− + 
2
pi + 
2
+ = 1.
Calculating the values for CAB,ˆ
√
(2J ′ + 1) is straight forward and included
in a summary of all single photon Rabi frequencies depending on polarization from
the states |0〉 and |1〉 in the Table 3.2.2.2. The table takes advantage of the fact
that the D1 manifold (2S1/2 → 2P1/2) and D2 manifold (2S1/2 → 2P3/2) have the
property that Γ0 ≡ γD1
√
I
2IsatD1
≈ γD2
√
I
2IsatD2
, differing only by about 3%. This
approximation will be used throughout to simplify expressions.
56
ΓAB for the D1 Manifold (
2S1/2 → 2P1/2)
F=0 F=1
mF 0 -1 0 1
A
|0〉 0 − 1√3Γ0 pi 1√3Γ0 + 1√3Γ0
|1〉 pi 1√3Γ0 − 1√3Γ0 0 +−1√3Γ0
ΓAB for the D2 Manifold (
2S1/2 → 2P3/2)
F=1 F=2
mF -1 0 1 -1 0 1
A
|0〉 −
√
2
3
Γ0 pi
√
2
3
Γ0 +
√
2
3
Γ0 0 0 0
|1〉 −−1√6Γ0 0 + 1√6Γ0 − 1√2Γ0 pi
√
2
3
Γ0 +
1√
2
Γ0
Table 3.1: Table of single photon Rabi frequencies. A summary of all single photon
Rabi frequencies, ΓAB, from the states A = |0〉 and A = |1〉 to both the D1 manifold
and the D2 manifold depending on polarization. The excited states are denoted by
their total angular momentum, F, and its projection on the z-axis, mF .
57
◦
320 330 340 350 360 370 380
2.×10-6
4.×10-6
6.×10-6
8.×10-6
1.×10-5
Laser Wavelength (nm)P
ro
b.
of
S
po
nt
an
eo
us
E
m
is
si
on
2P3/2 2P1/2
3[3/2]3/2
Figure 3.7: Probability of a spontaneous emission during a pi pulse as a function of
wavelength. Important to note is that missing factors of 2 have been corrected from
previous calculations, making it 4 times smaller than predicted prior.
If we only consider the D1 and D2 lines and ignore the hyperfine structure
since it only modifies the answer on the order of 2ωHF/∆, which in the worst case is
7× 10−4, then the total spontaneous emission rates for {|0〉 , |1〉}, reduce to simply
γ(p,|0〉) = γ(p,|1〉) =
Γ20
4
[
1
3
γD1
∆2
+
2
3
γD2
(∆′)2
]
(3.8)
where ∆ is the detuning from the D1 manifold and ∆′ is the detuning from the D2
manifold.
There is actually an additional energy level that must be considered for this
calculation due to its proximity, namely the 3[3/2]3/2 state, which is only 16.9 THz
away [69]. A good explanation of what a bracket state is and how to understand it
has been done by Senko [50] and won’t be repeated here. Writing the state in the
58
LS basis,
∣∣3[3/2]3/2〉 = √10
21
∣∣2P3/2〉+√ 8
21
∣∣4P3/2〉+ 1√
35
∣∣2D3/2〉− 2√
35
∣∣4D3/2〉 (3.9)
we see that the overlap with
∣∣2P3/2〉 is almost half. The lifetime of this state has
not been measured but has two separate calculated lifetimes of τ = 120 ns [70]
and τ = 150 ns [71] which gives a line-width of γ/2pi = 1.1 ± 0.3 MHz. If we
incorporate this bracket state as if it were simply, to leading order, a 2P3/2 state,
then the spontaneous emission in Eq. 3.8 is modified to
γ(p,|{1,0}〉) =
Γ20
4
[
1
3
γD1
∆2
+
2
3
γD2
(∆′)2
]
+
Γ2Bra
4
2
3
γBra
(∆Bra)2
(3.10)
where ΓBra is the single ion Rabi frequency for the bracket state and ∆Bra is the
laser detuning from the state. Using this equation, we can examine the probability
of spontaneous emission during a pi pulse as a function of wavelength which is shown
in Fig. 3.7. By measuring the Rabi frequency at the ions, we can also calculate what
the spontaneous emission probability is as a function of time. For a single Raman
beam and one ion, even out at 50 ms, there is still less than 4% probability. Still,
since each ion has an independent probability, p, to spontaneously emit, then the
probability to not have a single spontaneous emission event in a chain of N ions is
(1 − p)N . Using the same parameters, the probability of a single ion in a chain of
N ions spontaneously emitting is shown for multiple N in Fig. 3.8. Assuming a 10
ms total experimental time, Fig. 3.8 exhibits what this exponential scaling implies:
59
0 10 20 30 40 50
0.0
0.1
0.2
0.3
0.4
0.5
Time (ms)
P
ro
ba
bi
lit
y
of
a
S
in
gl
e
E
m
is
si
on
1 Ion
10 Ions
30 Ions
100 Ions
Figure 3.8: Probability of a single spontaneous emission in N ions. Due to the
exponential scaling with system size, spontaneous emission errors become a signif-
icant consideration at the time-scale of experiments at 30 ions. With 100 ions,
spontaneous emission will be a dominant error-source.
while for 10 spins and maybe even 30 spins this spontaneous emission can be treated
perturbatively, at 100 spins one of the dominant errors will be spontaneous emission.
We can use the exponential scaling of the spontaneous emission probability to
experimentally measure the spontaneous emission rate. With a system of ten ions,
we perform an experiment similar to that of Fig. 3.5a, except this time, instead of
letting the ions sit in the dark, we instead apply a single Raman beam to the chain
and wait. Any ions that measure as bright beyond the normal background have
had a spontaneous emission event. Using the camera, we can count the number
of ions bright and measure the rate at which the zero bright population decreases.
Just such a measurement is shown in Fig. 3.9 for ten ions. A baseline with no
applied light was measured and showed no change for the duration of the experiment.
Whereas when the light is applied, by fitting the decay of the 0 bright population
to an exponential of e−Nt/τ , we can extract the spontaneous emission rate, namely
γ(p,|{1,0}〉) = 1/τ = 2pi103 ± 6 mHz. This measured rate is consistent with the
60
● ● ● ● ● ● ● ● ● ● ● ● ●
■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■
■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆
0 5 10 15 20 25 30
0.0
0.2
0.4
0.6
0.8
1.0
Time (ms)
P
ro
ba
bi
lit
y ● 0 Bright■ 1 Bright◆ 2 Bright
Fit
Figure 3.9: Measured spontaneous emission rate. Experiment is to prepare N = 10
ions dark and apply a Raman beam for some time t and then measure. By using
the camera, we count the number of ions bright at the end of the experiment. We
fit the decay of the zero ions bright to an exponential, e−Nt/τ . The spontaneous
emission rate, γ(p,|{1,0}〉), is then calculated γ(p,|{1,0}〉) = 1/τ = 2pi103± 6 mHz
predicted value from Eq. 3.10 and Fig. 3.7.
Further, the spontaneous emission is difficult to improve beyond the observed
value. From Fig. 3.7, it is clear 355 nm is a local minimum in spontaneous emission
and in fact does not improve as the wavelength gets larger than 369 nm due to the
large drop in Rabi frequency. The spin-spin interaction used to generate entangle-
ment is proportional to the square of the intensity I, namely J ∝ (ηΩ)2/δ where
Ω ∝ I and δ is the detuning from the motional mode. Initially, this indicates that
increasing the intensity would effectively suppress the spontaneous emission, since
γ(p,|{1,0}〉) ∝ I and J is increasing quadratically with intensity, thus performing the
simulation faster than the spontaneous emission events. Yet in the slow gate regime,
where δ must be much greater than ηΩ, the optimal detuning is δ ≈ 3ηΩ. Plugging
this back in to the equation for J , we see that J is effectively proportional to the
intensity and that there is no suppression of the spontaneous emission in the slow
61
gate regime. As we scale up to larger systems, we will have to be more aware of
spontaneous emission in the simulations, especially since there is a 2/3 probability
that the event will take the atom outside the qubit basis.
3.3 State Detection
As discussed above, state detection is performed using state-dependent fluo-
rescence which is collected on either a photo-multiplier tube (PMT) or an intensi-
fied charge-coupled device (ICCD). The collection optic for the light is a 0.23 NA
microscope objective (CVI UVO-20.0-10.0) which is used in a finite conjugate con-
figuration with an effective front focal length of ∼ 18 mm, which is then imaged
by a doublet onto the measurement device, giving an effective total magnification
of ∼ 200. Between the objective and the camera is an adjustable zero aperture
iris (ThorLabs SM1D12CZ) for spatial filtering and two interference filters for color
selectivity at 370 nm, one which is a 10 nm narrow band around 370 nm (Semrock
FF01-370/10-25) and the other which is a long pass filter (Semrock LP02-355RS-
25). These filters effectively suppress the background scatter from room light to less
than 1 photon at the camera. We can switch from the camera to a PMT by using
a simple flipper mirror, where there is a an additional 370 narrow band filter for
better 355 nm suppression at the PMT. Since the PMT area is smaller than that of
the camera, there is a demagnifying lens in the PMT arm of the imaging system.
The PMT (Hamamatsu H10682-210) has a quantum efficiency of ∼ 30% at
370 nm. We use the PMT for diagnostics and calibrations because of the near
62
instantaneous readout and low noise. Typical mean photon counts on the PMT
within an 800 µs window is ∼ 18 while background counts are ∼ 3. The excess
background counts are due to reflections of the resonant light off of the Raman
optics after the trap reflecting back into the imaging optics. While not ideal, since
all data is taken using the ICCD camera which has spatial discrimination, the large
background is not problematic.
The ICCD camera (Princeton Instruments PI-Max3:1024i) has a smaller quan-
tum efficiency of ∼ 20% at 370 nm. While this is comparable to the PMT, there
are additional complications in any amplified CCD that decrease the signal-to-noise
ratio. Since the ICCD amplifies the photo-electrons by using an electron cascade,
the pixel values end up having a super-Poissonian distribution which has twice the
variance of the photon distribution [72]. This effectively halves the signal-to-noise of
the camera in single shot measurements due to a much larger overlap with the dark
distribution. In order to recover a comparable signal-to-noise, the camera exposure
is doubled to 1.5 ms. The drawback is that due to the off-resonant coupling of the
detection light to the F = 1 states of the 2P1/2 manifold, there is a greater chance
of misdiagnosis of dark as bright.
The PI-Max3:1024i has also presented significant technical difficulties. The
largest problem is that the charge transfer operations of the CCD are not efficient
when the camera run with a fast clock-speed at low light levels. This leads to a
displacement and elongation of the signal in the vertical shift direction. In order to
collect the signal, we have to integrate most of the CCD in the vertical direction
(approximately half). By using the hardware binning feature, which for a binning
63
&KDSWHU *HQHUDO2SHUDWLRQ)DFWRUV 
Full Frame Binning  
)LJXUHVKRZVDQH[DPSOHRIELQQLQJZLWKGXDOSRUWRSHUDWLRQIRUDQLQWHUOLQH
DUUD\HDFKSL[HORIWKHLPDJHGLVSOD\HGE\WKHVRIWZDUHUHSUHVHQWVSL[HOVRIWKHDUUD\
)LJXUHVKRZVWKH/LJKW)LHOGRegion of InterestH[SDQGHUVHWXSIRUIXOOVHQVRU
ELQQLQJ)LJXUHVKRZVDQH[DPSOHRIVLQJOHSRUWRSHUDWLRQIRUDIXOOIUDPHDUUD\

Figure 17. Dual Port Readout: 2 × 2 Binning of Interline CCD  

Figure 18.  Dual Port Readout: LightField Settings for 2 × 2 Binning of Interline CCD  
Figure 3.10: Schematic representation of CCD hardware binning. When using the
hardware binning on a PI-Max3:1024i, the photo-sensitive pixels are first shifted to
the adjacent storage pixels. For a binning of 2× 2, first two successive vertical shift
operations are performed, followed by two horizontal shift operations in the readout
register. Then once the four pixel’s worth of charge has been combined and shifted
to the output node, the charge is readout, creating only one “super-pixel” of data.
This is then repeated. Compliments of Princeton Instruments.
of m× n collects m pixels horizontally and n pixels vertically prior to read out, we
can generate “super-pixels” which recover most of the signal. The optimal hardware
binning parameters is a binning of 4 × 64, which is highly asymmetric in order to
compensate for the poor fidelity of the shift operations.
Another difficulty stemming from the hardware is that the base background
value of each pixel was dependent on the time between readouts. What we observed
is that for long experiments (∼ 5 ms or more), the background values of each pixel
grows to be comparable to the original signal size. The growth of this background
is due to the accumulation of dark charge which is exacerbated as the camera heats
up during heavy data collection. Thankfully, the photon signal is additional to this
background value and not competing with it, meaning that the effective signal to
64
noise is approximately the same, just with a larger offset. This issue made any
discrimination of the ion signal based directly on the pixel values prone to errors.
One final issue is the overlap of the point-spread functions (PSF) of the ion
fluorescence on the camera. With only an NA=0.23, the final width of the PSF when
observed on the camera is comparable to the effective inter-ion spacing, meaning
that ion crosstalk is a non-negligible source of error. This issue is compounded as
the axial trap is tightened or more atoms are added to the trap.
We examined multiple methods of analysis [73] and due to these issues, we
settled on a method of analysis that attempted to mitigate the effect of most of the
issues.
3.3.1 Camera Analysis
In order to reduce data handling overhead, the pictures are summed down to
a single row, reducing the dimension of the data to 1-D. By calibrating the center
of each ion PSF and its width on the camera, a sum of N Gaussians is fit to the
data where there are only N free parameters, namely the height of each Gaussian
distribution. A Gaussian is used as the fundamental distribution because the PSF
of a point source, namely an Airy disk, has a large overlap with a Gaussian. A
sum of Gaussians is used so that the ion crosstalk is diminished since the overlap of
distributions is taken into account by the fitting itself. There is an additional free
parameter which is the global offset to the whole distribution. This global offset
compensates for the growth of the background and only requires that there is at
65
least an ion width of dark area outside of the ion chain on each side. Once the data
has been fitted, the state of each ion reduces to the discrimination of the Gaussian
amplitude.
The fitting of this distribution requires prior knowledge of the ion number, po-
sitions and widths. Discrimination requires high statistics in order to determine the
optimal threshold and what the resulting errors are. We calibrate these parameters
prior to data collection by using 1000 images of an all bright chain and 1000 images
of an all dark chain. These calibration images are first averaged with the averaged
dark image being used as a background. The resulting background-subtracted dis-
tribution is used to determine the parameters for each ion’s Gaussian distribution.
An example of this averaged distribution can be seen in Fig. 3.11a. Once the Gaus-
sian parameters are calibrated, the Gaussian distribution is then fit to each of the
2000 calibration images. A histogram of each set of images is generated in order to
optimize the threshold.
In the ideal case, the point where the bright histogram and dark histogram
have equal error is the best point to discriminate. In the current apparatus, this is
not the case because despite the fact that the fitted Gaussians take into account the
PSF overlap, the photon counts are low enough that the data is fairly stochastic.
Figure 3.11b shows a single image of bright 24 ions which has only been background
subtracted and shows how large signal spikes can occur even in-between ions. This
results in a residual inter-ion crosstalk which increases the bright and dark histogram
separation biasing the discrimination towards the bright state. This is especially
apparent when looking at spin chain configurations with alternating spin states. To
66
0 50 100 150 200 250
0
500
1000
1500
2000
Column Number
Sig
na
lH
eig
ht
(a.u.)
Data
Fit
Thresholds
0 50 100 150 200 250
0
200
400
600
800
Column Number
S
ig
na
lH
ei
gh
t(a.u.
)
Avg Signal
Thresholds
(a)
(b)
Figure 3.11: Example of camera summed row data. Averaged signal from a twenty-
four ion calibration (a). Ions are optically pumped dark for 1000 images, which are
averaged together to form the background. Ions are rotated bright for 1000 images,
which are then averaged and background subtracted, producing an average signal.
This is used to calibrate fit of the distribution. Also displayed are the thresholds
associated with each ion, indicating height that fitted Gaussian must be in order to
be determined bright. Single image summed row distribution (b). Displayed here
are the background subtracted single image data and the associated Gaussian fit.
Small photon numbers lead pixel heights to have a wide variance. Also shown are
threshold Gaussian amplitudes, showing that threshold sizes are reasonable.
67
1 5 10 15 20 24
0.87
0.88
0.89
0.90
0.91
0.92
0.93
0.94
Ion Num
B
rig
ht
S
ta
te
Fi
de
lit
y Average Fidelity: 92%
1 5 10 15 20 24
0.975
0.980
0.985
0.990
Ion Num
D
ar
k
S
ta
te
Fi
de
lit
y
Average Fidelity: 98.5%
(a)
(b)
Figure 3.12: Example raw ion fidelities using ICCD. Raw dark/bright state fidelities
after discrimination (a)/(b). Due to bias required in bright state thresholds to
suppress inter-ion crosstalk, the resulting raw bright state fidelities are typically
worse than raw dark state fidelities. For this string of twenty-four ions, the averaged
dark state fidelity is 98.5% and the averaged bright state fidelity is 92%
68
compensate for this, we apply thresholds to the data with the following equation:
T (I) = T0 + ηB(I) (3.11)
where T is the threshold value for ion I, T0 is a threshold offset chosen to be slightly
larger than the individual summed column noise, and η is a scaling coefficient for the
background due to adjacent ion overlap B(I). The calculated threshold is show in
Fig. 3.11a. In comparison to the relative background, it is clear why this threshold
is scaled by background from the other ions. This type of operation allows for the
maximal sensitivity to staggered states. The drawback is that this greatly increases
the error on the bright state. In Fig. 3.12, we show the resulting bright and dark
fidelities, where there is an average of fidelity of 98.5% for the dark and 92% for the
bright.
Ultimately these fidelities are insufficient to make any observations of any
global spin observable as the observable fidelity F ∼ (f)N where f is the single ion
fidelity in an N spin system. For a twenty-four spin system with a single ion fidelity
of 95%, the global observable fidelity is only 29%, where now the single largest error
is just due to detection error.
3.3.2 Detection Correction
The solution to this difficulty is to leverage the knowledge we have of the
detection infidelities in order compensate for the detection error. These detection
errors are corrected using a data processing technique described in [74]. We will
69
first lay out an intuitive understanding of this error correction and then describe
the technique in [74] as it is implemented in our analysis.
For a single qubit, the most generic representation is a 2×2 matrix:
D =
[
1− p0 p1
p0 1− p1
]
(3.12)
where p0 (p1) is the error on the bright (dark) state. If we take some ideal input
probability distribution of gi, then the observed probabilities fi are F = DG. This
is a just a prediction of what the observation will be based on our knowledge of p0
and p1. In order to recover the ideal probabilities from our imperfect measurement,
we just need to invert the equation G = D−1F . This works very well for a single
ion.
The key observation of [74] was that for an N spin system, the 2N×2N matrix
Mji equivalent to the single spin matrix D, where fj =
∑2N
i=1 Mjigi, has a very simple
tensor product structure. In other words, we can instead express Mji as
M = ⊗Nk=1Dk (3.13)
where Dk is the single ion detection matrix for ion k. This is critical, as a general
matrix M is not necessarily invertible and even if it is, M has 22N elements which
makes constructing it difficult, let alone inverting it due to the cost in computa-
tional resources as N gets larger. Since M can be expressed as this tensor product,
inverting M is as simple as inverting each Dk which has a very simple analytic
70
expression:
M−1 = ⊗Nk=1D−1k = ⊗Nk=1
[
1− p′0,k p′1,k
p′0,k 1− p′1,k
]
(3.14)
where
p′0,k =
p0,k
p0,k + p1,k − 1
p′1,k =
p1,k
p0,k + p1,k − 1
(3.15)
with p0,k (p1,k) is the bright (dark) error for ion k. Solving for the real probability dis-
tribution gi from the observed probability distribution fi, we find gi =
∑2N
j=1 M
−1
ij fj,
allowing us to recover the real probabilities by just constructing M−1.
Now in general M−1 is still a very large matrix and computationally costly
to compute. For some idea of the resources involved, a twenty-four spin M−1 has
248 ∼ 2.8 × 1014 elements and is not sparse. If each entry is an 8-byte float, the
amount of memory required to just store the matrix is 2.2 petabytes.
One additional observation in [74] greatly improves this resource cost. Assume
that for an N spin system, we want to measure some observable Oˆ which is some
combination of spin projections. We can write Oˆ as a tensor product of Pauli
matrices, namely Oˆ = ⊗Nk=1σµkk , where σµkk is a component of the Pauli matrices
for spin k with µk = 1, 2, 3 ≡ x, y, z and µk = 0 is the identity matrix. In order to
measure a particular axis of the Bloch sphere, that axis is coherently rotated into the
zˆ basis which is the measurement basis. Thus Oˆ can be written as diagonal matrix
with Oˆ = ⊗Nk=1diag (σµkk ) where diag (σµkk ) = [1,−1] for µk 6= 0 and [1, 1] for µk = 0,
implying that Oˆ can be simplified to a vector of length 2N . With a probability
71
distribution of gi, the expectation value of Oˆ becomes simply 〈Oˆ〉 =
∑
i Oˆigi. Yet,
applying the results above, we find
〈Oˆ〉 =
∑
i
Oˆigi
=
∑
i
Oˆi
∑
j
M−1ij fj
=
∑
j
(∑
i
OˆiM
−1
ij
)
fj
=
∑
j
Oˆcjfj
(3.16)
where Oˆc now defines a corrected operator. This corrected operator can be written
using M−1 = ⊗Nk=1D−1k , namely Oˆc = ⊗Nk=1[diag (σµkk )D−1k ]. This equation reduces
the computation of a 2N × 2N matrix to only computing a 2N vector. The elements
for this vector are expressed simply as
Oˆc = ⊗Nk=1
[
(1− 2p′0,k),−(1− 2p′1,k)
]
. (3.17)
Generating this vector for Oˆc can be accomplished very efficiently, making this an
effective technique even at larger ion numbers (N < 30).
We implement this in the data analysis by first taking the discriminated spin
values and constructing a probability vector of all 2N possible spin states for the
chain. Then the ion errors from the camera calibration are used to either generate
the full matrix M−1 when we are measuring the global state probabilities or only the
single spin magnetization vector Oˆck of the kth spin for magnetization measurements.
72
These single spin magnetization vectors can be multiplied together to generate cor-
relation measurements to as many orders as required. The distinction of purpose
arises since when constructing the corrected operators, the observable is being cor-
rected, not the probabilities. When using the full matrix M−1 on the measured
probability vector, the matrix effectively moves population around to compensate
for the observation error whereas the corrected observable simply rescales the final
value of the observable. Ultimately, constructing the full matrix and finding the
real probabilities gi and then constructing the observable versus constructing the
corrected observable result in slightly different answers because the full matrix is
more sensitive to the detection error calibration.
If we let the bright and dark detection errors be the same p0 = p1 = p for
the sake of brevity, then a relative error e in the calibration of p is δp/p ∼ e. The
support of Oˆ, np, is equivalent to the number of Pauli matrices in the tensor product
of Oˆ. Thus the single spin magnetization has a support of np = 1 whereas a two-
point correlation function has a support of np = 2. Now, the error in the corrected
observable Oˆc due to calibration error is δ〈Oˆc〉/〈Oˆc〉 ∼ 2nppe. Constructing the
matrix has a support of np = N , whereas constructing Oˆ
c only has the support of
the observable. For this reason and also since it is computationally more tractable,
it is always preferable to use the corrected observable.
73
3.3.3 Independent Component Analysis (ICA)
As our working system sizes have increased, there has been a dramatic increase
in the computation time required to perform the data analysis. This slow down is a
direct consequence of the increase in the number of free parameters in the Gaussian
amplitudes fit. For a system of 26 spins, analysis of a scan with 21 points and
1000 experiment per point takes over 24 hours. This is not tenable, especially when
considering system sizes of > 50 spins and the corresponding non-linear increase in
analysis time. Any replacement of this Gaussian amplitude fit needs to not only be
more computationally efficient, but also automatically determine the ion positions
on the camera.
The ideal case would be to generate a set of basis images where only one ion is
bright and the rest are dark. Analysis would then just be the convolution of the basis
images with the data images, where the basis image would act as a mask, selecting
out the fluorescence from only a one ion, leading to a signal value for that ion. This
would then be discriminated in the exact same way as the Gaussian amplitudes.
Further, the basis images could be used to extract and compensate the inter-ion
crosstalk, further improving over the Gaussian amplitude difficulties described.
The difficulty is that generating these basis images is arduous as preparing a
single spin excitation for each spin has a high experimental time cost, whether it be
using single-site addressing or by moving a single ion around. A better method would
be to use an algorithm to identify what fluorescence belongs to each ion. One very
good possibility is to use independent component analysis (ICA). ICA is an ideal
74
Figure 3.13: Fast ICA blind source separation. Fast ICA is used here to identify
super-imposed sources in noisy data. First image is the “observed” signals. Second
image is the true sources used to generate the data. Third are the components
recovered by the ICA algorithm. These recovered signals are the same as the true
sources up to a minus sign. All credit for image to scikit-learn.
algorithm for separating superimposed signals, one classic example being the “blind
source separation” problem. In the blind source separation problem, there are, for
example, three audio sources which are recorded by three separate microphones.
The problem is to recover the true sources by examining their superposition from
the three separate recordings, which ICA does remarkably well.
We can map the problem of identifying ions directly to the blind source sep-
aration [75]. This is done by reducing each image to a 1-D array where each row
is appended in sequence. Finding the ions in the image now is the same as find-
ing its unique “frequency” in the waveform of the picture. As long as the number
of images input is much greater than the number of ions, this now reduces to the
same problem as that of blind source separation. Now there is one very stringent
75
Figure 3.14: Basis images from FastICA algorithm. “FastICA” from the is applied
to 10,000 images of ions with pi/2 rotation in order to make signals uncorrelated.
Result is 10 images, each with a single ion separated out. These basis images are
then used as a mask for the data to determine ion brightness.
requirement to using ICA that can make or break the algorithm and that is the
sources being examined must be statistically independent. For the ion system, this
requires essentially that 〈σizσjz〉 = 0 for all pairs. This is achieved by simply rotating
all spins by pi/2. An example of this being applied to a 10 ion system can be seen
in Fig. 3.14. The particular implementation of ICA is from the open-source Python
package of Scikit-Learn and called “FastICA”.
While this algorithm works well for small system sizes (∼ 10 ions), when the
system gets much larger the algorithm has a high failure rate. Initially, the quality
of the images were suspected, but later it was realized that in fact the pi/2 pulses on
the outer ions were no longer pi/2 but were under-rotated. Thus their signals were
no longer statistically independent and so the ICA could no longer appropriately
distinguish them as a source. In order to implement this as a replacement to the
current analysis, it requires that the pi/2 pulses are even across the chain. Currently
76
the Raman beams are used to perform the pi/2 pulse leading to the issues described,
but this can be circumvented by using microwaves for the calibration images which
should produce much more even rotations.
77
Chapter 4: Fourth Order Stark Shift
The pervasive challenge facing all quantum information platforms is the unde-
sired interaction of the qubit with environment. In trapped ions, one such coupling
to the environment is via the modulation of the qubit energy splitting by stray
magnetic fields. This can be circumvented by using clock-states, allowing for co-
herence times exceeding 10 minutes [76, 77]. However, the use of clock-state qubits
by definition does not easily allow the direct generation of certain classes of Hamil-
tonians that are equivalent to the modulation of qubit energy splittings [58]. In
quantum computing, such control is desirable for efficiently realizing universal logic
gate families such as arbitrary rotations [78].
In this chapter, I propose and demonstrate the use of a fourth-order Stark shift
to achieve fast, individually addressed, single-qubit rotations in a chain of 171Yb+
ions. We experimentally realize up to a 10 MHz shift on the qubit splitting with
only moderate amounts of laser power. We exploit this control in a quantum system
of 10 ions by preparing arbitrary initial product states and applying an independent
programmable disordered splitting on each lattice site in a quantum simulation, all
demonstrated with low cross-talk.
78
4.1 Fourth-order Stark Shift Theory
The studies reported here are performed on a linear chain of 171Yb+ ions,
but can be generalized to any species of clock qubits. Here we will adopt the
notation 2S1/2 |F = 0,mf = 0〉 and 2S1/2 |F = 1,mf = 0〉 are denoted as |0, 0〉 and
|1, 0〉 respectively.
We assume the ions are irradiated using an optical frequency comb generated
from a mode-locked laser with a center frequency detuned by ∆ from the 2P1/2
manifold and by ωF − ∆ from the 2P3/2 manifold. We assume that the pulse area
of each laser pulse is small and has only a modest effect on the atom, and that
the intensity profile for each pulse is well approximated by a hyperbolic secant
envelope [64]. Under these assumptions, the kth comb tooth at frequency kνrep
from the optical carrier has a resonant S → P Rabi frequency [79],
gk = g0
√
piνrepτsech(2pikνrepτ) (4.1)
where τ is laser pulse duration, g20 = γ
2I¯/2I0, I¯ is the time-averaged intensity of the
laser pulses, I0 is the saturation intensity of the transition, and γ is the spontaneous
decay rate. Since
∑∞
k=−∞ g
2
k = g
2
0, and assuming the parameters specified above,
the second-order Stark shift E
(2)
α of state |α〉 due to the frequency comb can be
79
computed for an arbitrary polarization (taking ~ = 1) [64,66]:
E
(2)
00 =
g20
12
(
1
∆
− 2
ωF −∆
)
E
(2)
10 =
g20
12
(
1
∆ + ωHF
− 2
ωF − (∆ + ωHF )
)
.
(4.2)
Here we neglect all excited state hyperfine splittings since they only contribute to
the Stark shifts at a fractional level of ∼ 10−5. We also ignore all other states
outside of the P manifold since their separation from the ground S states are too
far detuned from the applied laser fields to give appreciable Stark shifts.
Assuming that 20 mW of time-averaged power is focused down to a 3 µm waist,
the differential second-order Stark shift on the qubit splitting is δω(2) = E
(2)
10 −E(2)00 =
−1.5 kHz.
We will show that there is a fourth-order effect that can be much larger than
the differential second-order Stark shift when using a frequency comb for specific
polarizations of the beam. An intuitive understanding can be gained by considering
that any two pair of comb teeth, k0 and k1, have a beat-note frequency (k0 −
k1)2piνrep. If the bandwidth of the pulse is large enough, then there will be beat-
notes that are close to the ground state hyperfine splitting. Assuming that none
are on resonance, these off-resonant couplings can have a large effect on the ground
states, as much as three orders of magnitude larger than the differential AC Stark
shift.
We first calculate the fourth-order Stark shift in the simplified case of just two
comb teeth and one excited state of the 171Yb+ level structure (see Fig. 4.1), equiv-
80
1,1
0,0
1,0
1,-1
∆
ωHF2S1/2
2P1/2
2P3/2
ωZee
ωF
Γ0,ω0
Γ1,ω1
δ
Figure 4.1: Schematic representation of the electron energy levels of 171Yb+. When
two phase-coherent colors of light are applied to the atom which have a beatnote
approximately equal to the qubit splitting, there is an effective fourth-order differ-
ential light shift which can be much larger than the second-order differential Stark
shift.
81
alent to two phase coherent continuous wave (CW) beams in a three level system.
Let the excited state |e〉 have frequency splitting ωe from the |0, 0〉 ground state,
and the absolute frequencies of the comb teeth k0 and k1 be ω0 and ω1 respectively.
Also, let the polarization of each tooth, i, be defined as ˆi = ˆ = −σˆ− + 0pˆi+ +σˆ+
with |−|2 + |0|2 + |+|2 = 1 where σˆ−, pˆi, and σˆ+ are the polarization basis in the
frame of the atom. In the rotating frame of the electro-magnetic fields of the laser,
we can write the Hamiltonian
H =H0 + V
=δ |1, 0〉 〈1, 0|+ ∆ |e〉 〈e|
+
Γ0
2
|0, 0〉 〈e|+ Γ
1
2
|1, 0〉 〈e|+ h.c.
(4.3)
where H0 contains the diagonal terms and V includes the off-diagonal terms induced
by the laser, δ = ωHF − (ω0−ω1), Γi = g0C(ˆi) is the resonant Rabi frequency from
beam i with a dipole coupling matrix element C(ˆi) for polarization ˆi. The fourth-
order correction E
(4)
n to the ground state energy levels, from perturbation theory,
has the following form:
E(4)n =
∑
j,l,m6=n
Vn,mVm,lVl,jVj,n
En,mEn,lEn,j
− |Vn,j|
2
En,j
|Vn,m|2
(En,m)2
− 2Vn,nVn,mVm,lVl,n
(En,l)2En,m
+ V 2n,n
|Vn,m|2
(En,m)3
.
(4.4)
Here j, l,m, and n each represent different energy levels, Va,b = 〈a|V |b〉, Ea,b =
E
(0)
a − E(0)b is the unperturbed energy difference between the states |a〉 and |b〉.
Applying this to the Hamiltonian above, the last two terms are zero since V has no
82
Galilean Telescope:
Mag = 3x
AOD Imaging System
ICCD Camera
369nm 
355nm 
Ion Imaging Objective System
Dichroic Optic
NA 0.23
Figure 4.2: Diagram of optics imaging 355nm light onto single ion. Final spot size
is < 3 µm, giving rise to controllable and individual-addressed Stark shifts on the
qubits. This optical system utilizes a NA 0.23 objective lens for state detection of
the ions at 369nm. Since the AOD is not imaged, deflections at the AOD correspond
to displacement at the ions. This maps RF drive frequency to ion position, enabling
control of the horizontal position of the beam.
diagonal terms leaving the fourth-order Stark shifts of the qubit levels,
E
(4)
00 =−
|Ω|2
4δ
E
(4)
10 =
|Ω|2
4δ
.
(4.5)
In these expressions, we assume δ  ∆ and Γ0 ∼ Γ1. We also parametrize Ω =
Γ0Γ1/2∆, which is the resonant (δ = 0) stimulated Raman Rabi frequency.
The above derivation is valid for any three level system. We now include the
more complete case in 171Yb+ where all excited states with major contributions,
namely the 2P1/2 and
2P3/2 manifolds, are considered. Calculating the fourth-order
Stark shift on any state |n〉 reduces to computing its shift due to all other states cou-
pled via a two-photon Raman process by fields at frequencies ω0 and ω1. In
171Yb+,
this means we must consider all hyperfine ground states. The two Zeeman states,
|F = 1,mf = ±1〉, of the ground state manifold will be denoted as {|1, 1〉 , |1, -1〉}.
83
To calculate the fourth-order Stark shift, we sum over all states |a〉 6= |n〉,
E(4)n =
∑
a6=n
Ω2n,a
4δn,a
(4.6)
where Ωn,a is the two-photon Rabi frequency between |n〉 and |a〉, δn,a = ωa− (ω0−
ω1), and ωa = E
(0)
a −E(0)n . Computing all of the relevant Rabi frequencies Ωn,a under
the same assumptions as in Eq. 4.2 [66], we find
Ω00,10 =
(
0−
1∗
− − 0+1∗+
)
Ω0
Ω00,1-1 = −
(
0−
1∗
pi + 
0
pi
1∗
+
)
Ω0
Ω00,11 =
(
0+
1∗
pi + 
0
pi
1∗
−
)
Ω0
Ω10,1-1 =
(
0−
1∗
pi + 
0
pi
1∗
+
)
Ω0
Ω10,11 =
(
0+
1∗
pi + 
0
pi
1∗
−
)
Ω0.
(4.7)
Here Ω0 =
g20
6
(
1
∆
+ 1
ωF−∆
)
and g20 = γ
2I¯/2I0. From Eq. 4.7, we see that if ˆ =
σˆ±, the Rabi frequency Ω00,10 is maximized and equal to Ω0. If instead ˆ = βˆ ≡
1/2σˆ− + 1/
√
2pˆi + 1/2σˆ+, which corresponds to a circularly polarized input beam,
then Ω00,10 = 0 while all other Rabi frequencies are equal to Ω0/
√
2. It should be
noted that a linear polarization from a single beam cannot drive Raman transitions
between any of the 171Yb+ hyperfine ground states. These polarizations are the two
which provide the largest Rabi frequencies, while all others have smaller effective
Rabi rates, so we dwell on these two cases. An important note is that in the case of
ˆ = βˆ, E
(4)
10 = 0 because the shifts from |1, 1〉 and |1, -1〉 are equal and cancel each
84
other.
We now compute the differential fourth-order Stark shift on the qubit states
|1, 0〉 and |0, 0〉, δω(4) = E(4)10 − E(4)00 ,
δω(4) =

Ω20
2δ00,10
when ˆ = σˆ±
Ω20
8
(
1
δ00,11
+ 1
δ00,1-1
)
when ˆ = βˆ.
(4.8)
Finally, we generalize to incorporate all possible pairs of comb teeth. The
two-photon Rabi frequency for any two comb teeth k0 and k1, where k1 − k0 = l is
Ωn = gk0gk0+l/2∆ ≈ Ω0sech(pilνrepτ) [79]. Let j be defined such that |ωa − 2pijνrep|
is minimized, assuming that it is nonzero. If we now plug this into Eq. 4.6 summing
over all comb teeth, we find
E(4)n =
∑
a6=n
Ω2n,a
4
∞∑
k=−∞
sech2((j + k)piνrepτ)
δn,a − k(2piνrep)
=
∑
a6=n
Cn,a
Ω2n,a
4δn,a
(4.9)
where δn,a = ωa − j(2piνrep), and
Cn,a =
∞∑
k=−∞
sech2((j + k)piνrepτ)
1− k(2piνrep)/δn,a . (4.10)
Because the denominator in Eq. 4.10 grows rapidly with k, only the closest few
beatnotes are important, and as long as 2piνrep  ωZee, then E(4)10 remains zero for
85
ˆ = βˆ. The differential fourth-order Stark shift then becomes
δω(4) =

C00,10 Ω
2
0
2δ00,10
when ˆ = σˆ±
Ω20
8
(
C00,11
δ00,11
+ C00,1-1
δ00,1-1
)
when ˆ = βˆ.
(4.11)
Assuming the same parameters as with the second-order Stark shift (20 mW of
time-averaged power focused down to a 3 µm waist) and with νrep = 120 MHz and
τ = 14 ps, we find that the fourth-order shift is
δω(4)/2pi =

247 kHz when ˆ = σˆ±
132 kHz when ˆ = βˆ.
(4.12)
This result is ∼ 100 times larger than the differential second-order Stark shift for
the same parameters. Comparing the fourth and second-order expressions, we find
that δω(4)/δω(2) ∝ g20/(ωHF δ), clearly defining the regime where the fourth-order
shift dominates. The second-order shift only becomes larger with a hundred-fold
reduction in the laser intensity, corresponding to an applied shift below 10 Hz.
Since the differential fourth-order shift can easily be made large as shown above, it
is a practical means to control a large number of qubits.
4.2 Experimental Setup
The laser used to generate the fourth-order Stark shift is the same mode-
locked, tripled, ND:YVO4 used for coherent manipulation of the qubit (see section
86
3.1.4). The optical access of our current vacuum chamber restricts the polarization
of the Stark shifting beam since the magnetic field is orthogonal to all viewports,
prohibiting the use of pure σ± light. However, as discussed earlier, the differential
fourth-order Stark shift has two possible polarizations with large shifts for a single
beam: the first is pure σ±, the second is ˆ = βˆ ≡ 1/2σˆ− + 1/
√
2pˆi + 1/2σˆ+. We use
the βˆ polarization which slightly reduces the maximum shifts applicable, but does
not require pure σ±.
The small spot size required to individually apply a shift to each qubit is
achieved by using the imaging objective designed for qubit state readout. Since the
cycling transition of 171Yb+ is 369 nm and the center wavelength of the modelocked
laser is 355 nm, we use a Semrock dichroic beam combiner (LP02-355RU-25) for
separating the 355 nm laser from the resonant light at 369 nm (Fig. 4.2). Guided by
simulations of the optical system in the commercial ray-tracing software, Zemax [80],
we focus the 355 nm light down to a less than 3 µm horizontal waist using an
objective lens with a 0.23 numerical aperture.
In order to address each ion in a chain of up to 10 sites, we use an acousto-
optical defelector (AOD, Brimrose CQD-225-150-355). Since the AOD is not imaged,
it maps the rf drive frequency to ion position and the rf power of that drive frequency
to the applied intensity. The rf control is implemented using an arbitrary waveform
generator (AWG, Agilent M8109A), because it allows precise, easy, and arbitrary
control while being easily reconfigurable. The differential fourth-order Stark shift
is a direct change in the energy splitting of the qubit, so unlike in stimulated Ra-
man processes, phase coherence does not require optical phase stability or even rf
87
t0
t0
t0
t0
t0
System Dynamics
Nt0 Nt0 Nt0 Nt0 Nt0
Ion 1
Ion 2
Ion 3
Ion N-1
Ion N
Figure 4.3: Sketch of a typical raster pulse sequence. When the light is evenly
distributed across N ions, the applied fourth-order stark shift diminishes by 1/N2
due to the quadratic dependence on intensity. We recover a linear dependence on
ion number by rastering the beam, or applying a large shift for a short time, t0
sequentially to the ions. As long as each pulse chapter of length Nt0 is much shorter
than the interaction time-scale, then the shift on each ion is then proportional to
1/N .
phase stability, but only depends on the integrated time-averaged intensity. Thus
phase-coherent control only requires timing resolution better than the period of the
differential fourth-order Stark shift, which is easily achieved with the AWG. The
AWG also allows the application of many frequencies to the AOD, which will Stark
shift multiple ions simultaneously. Additionally, the AWG gives arbitrary amplitude
control of each frequency, providing time-dependent amplitude modulation of the
four photon Stark shift.
Due to the quadratic dependence of the differential fourth-order Stark shift on
intensity, when we divide the optical power across N ions, each ion’s fourth-order
88
Stark shift is diminished by a factor N2,
δω(4)(ion) = max(δω(4))/N2. (4.13)
In order to recover a linear dependence, we “raster”, or rapidly sweep, the beam
position from site to site across the chain. If this rastering occurs much faster
than the dynamics of the system, then the effective fourth-order shift can be safely
time-averaged, yielding
δω(4)(ion) = max(δω(4))
mt0
T
(4.14)
where m is the number of raster cycles in the total elapsed time T and t0 is the
time the light is applied to each ion in a single cycle. In order for the raster to be
fast enough to justify averaging the Stark shift, the length of each raster cycle, Nt0,
must be small compared to the total elapsed time T = Nmt0. We assume that any
time where the beam is not hitting any of the ions during a raster cycle is small
compared compared to the raster cycle duration and can be neglected. Substituting
into Eq. 4.14,
δω(4)(ion) = max(δω(4))
1
N
(4.15)
which recovers a linear dependence on the system size. In Fig. 4.3, we show a
diagram of an example raster sequence. The limitation on this technique is how
small t0 can be made. In our case, t0 is limited by the rise time of the AOD, which
is approximately 50 ns, which is still fast compared to N/δω(4) and very fast when
89
compared to a mechanical deflector rise time.
4.3 Experimental Demonstration
Using Ramsey spectroscopy [81], we measure the total Stark shift on the
qubit splitting from the applied light. A quadratic dependence on the intensity
distinguishes the fourth-order Stark shift from the typical linear dependence of the
second-order AC Stark shift (Eq. 4.11). By measuring the total shift as a function
of applied time-averaged optical power, the data in Fig. 4.4a demonstrates that the
observed shift is consistent with the I¯2 dependence of the fourth-order Stark shift.
By translating the ion through the beam, we measure the horizontal beam
waist by fitting the resulting Stark shift to the square of a Gaussian distribution
(Fig. 4.4b):
δω(4)(∆x) = δω(4)(0)
(
e−2∆x
2/σ2
)2
. (4.16)
We measure the horizontal waist to be σ = 2.68 ± 0.03 µm. This small waist
allows for independent control of qubits. In Fig. 4.5a, we show how qubit 5 can be
driven in a ten ion system with only minimal crosstalk of approximately 2% on the
adjacent spins (ions 4 and 6). In this configuration, the ions are separated by 2.76
and 2.64 µm respectively. By increasing the distance between ions, we can decrease
the crosstalk on adjacent spins. For example, in a system of two spins separated
by 7 µm, we individually drive each ion with the cross-talk ≤ 2× 10−5 over a time
t = 30× 2pi/δω(4).
As indicated above, the rf drive frequency maps to position at the ion chain,
90
0.00
0.05
0.10
0.15
0.20
0.25
Li
gh
tS
hi
ft
(M
H
z)
-2 -1 0 1 2
-4
-2
0
2
4
Ion Position (μm)
R
es
id
ua
ls
0
2
4
6
8
10
12
Li
gh
tS
hi
ft
(M
H
z)
0 50 100 150 200
-4
-2
0
2
4
Applied Optical Power (mW)
R
es
id
ua
ls
(a)
(b)
Figure 4.4: Measured Stark shifts as a function of optical power and position. Mea-
sured fourth-order Stark shift as a function of optical power with fit residuals (a).
We fit the measured light shift as a function of optical power to Eq. 4.11 for ˆ = βˆ
taking into account an astigmatism of the imaging optics resulting in the vertical
waist being ∼ 1.5 times the horizontal and find very good agreement showing that
the light shift arises from the fourth-order Stark shift. Measurement of the beam
waist at the ion with fit residuals (b). By translating the ion through the beam with
a fixed applied optical power of 40 mW, we extract the horizontal optical waist at
the ion. We found this to be 2.68 µm.
91
0 5 10 15 20 25 300.0
0.2
0.4
0.6
0.8
1.0
Time (
s)
Pro
ba
bili
ty
|1,0
-40 -20 0 20 400.0
0.2
0.4
0.6
0.8
1.0 -15 -10 -5 0 5 10 15

F (MHz)
Pro
ba
bili
ty
|1,0

X (m)
Ion 1
Ion 2
Ion 3
Ion 4
Ion 5
Ion 6
Ion 7
Ion 8
Ion 9
Ion 10
(a)
(b)
Figure 4.5: Resolution of individual addressing beam. Observed crosstalk of beam
applied to one ion (a). By applying light to only ion 5 in a chain of 10, we measure
the crosstalk on the nearest neighbors, ion 4 and 6, to be only 2%, which is consistent
with our measured horizontal beam waist and the ion separation. Solid line is a fit
to an exponential decaying oscillation with decay parameter τ = 133µs, which is a
2% error per pi-pulse. Individual ion signal as the beam is swept over a chain of ten
ions (b). By scanning the AOD drive frequency for a fixed power and duration, we
map the fourth-order Stark shift as a function of drive frequency. This corresponds
to a displacement of beam position at the ion chain. The effective scanning range
of the AOD is approximately 30 µm.
92
while the small spot size allows for individual control of the ions. In Fig. 4.5b, we
show this mapping in a chain of ten ions by scanning the drive frequency of the AOD
while fixing the rf power and time. The difference in the applied fourth-order Stark
shift of each ion is due to the rf bandwidth of the AOD, since the diffraction efficiency
is lower at the extremes of the bandwidth. In the current optical setup, a change of
10 MHz to the drive frequency corresponds to a displacement of approximately 3.4
µm along the ion chain.
This control enables the preparation of arbitrary, high-fidelity product states
when the individual addressing beam is used in conjunction with global qubit oper-
ations from the Raman beams. In Fig. 4.6, we illustrate a pulse sequence used to
generate a product state. This method, effectively a Ramsey sequence, is used to
prepare a spatially-alternating spin state, which is the most difficult state to produce
since it is the most susceptible to crosstalk. We observe a fidelity of 87± 1% for the
desired state, which includes all state preparation and measurement (SPAM) errors.
This fidelity is consistent with a 2% error of the pi-pulses on five of the ions arising
from the intensity noise and the small inter-ion crosstalk of the individual addressing
beam, and some residual infidelity stemming from the off-resonant coupling of the
ground states and the ion detection error.
Two features of the experimental noise should be noted. The first stems from
the quadratic nature of the fourth-order light shift. This quadratic dependence on
the intensity doubles the fractional uncertainty in the light shift relative to small am-
plitude noise in the laser intensity. The second arises from the off-resonant coupling
of levels in the ground state manifold. This off-resonant coupling leads to effective
93
time
R (�/2) (�/2)x
All Rx
All
R (�)z1
R (�)z3
R (�)z5
R (�)z7
R (�)z9
Ion 1
Ion 2
Ion 3
Ion 4
Ion 5
Ion 6
Ion 7
Ion 8
Ion 9
Ion 10
Figure 4.6: Pulse sequence for preparing a string of 10 ions in a staggered spin
configuration. All 10 ions are prepared in |0, 0〉 and then a global pi/2 pulse is
applied. Depending on the state being prepared, some number of the ions have a
pi phase shift applied, creating the desired configuration. A final global pi/2 pulse
projects the configuration back into qubit basis, completing the effective Ramsey
sequence.
Rabi dynamics that results in a mixture of the qubit states and the coupled states
that is ∼ 1
2
Ω2n,a/(δ
2
n,a + Ω
2
n,a). In the experiment, light shifts of order 5 MHz are
therefore expected to mix qubit state and the coupled state up to 20%. Mixing the
states together is effectively dissipation, causing decoherence of the spin. For this
reason, during coherent operations we restrict the shift size to < 300 kHz, where
this mixture is less than 2.5% the qubit population.
The primary limitation in the current apparatus is the intensity applied to
each ion, especially those on the edge of the chain due to the bandwidth of the
AOD. The maximum intensity on each ion is simply
Iion =
2piP¯ (NA)2
λ2
DEν (4.17)
94
where P¯ is the time-averaged power into the AOD, NA is the objective numerical
aperture, λ is the wavelength of the light, DEν is the diffraction efficiency of the AOD
at the drive frequency, ν, corresponding to the ion position. By enlarging the NA
of the objective lens, the intensity applied to each ion would greatly increase while
simultaneously lowering the inter-ion crosstalk. Further, improving the diffraction
efficiency and bandwidth of the AOD will allow more ions to be addressed. By
implementing changes on both of these elements, we should be able to address 20+
ions without difficulty.
4.4 Applications of the Fourth-order Stark Shift
The freedom and control afforded by an individually addressed, Stark-shifting
beam opens many possibilities that were previously inaccessible to clock state qubits.
One such new application is that we can now apply site dependent transverse mag-
netic fields to an interacting Ising spin system [40]. Since the strength of each field
is controlled by the rf amplitude from the AWG, we are able to quickly generate
hundreds of different random instances of individual ion fields in a reproducible way.
Schematic representations are shown in Fig. 4.7.
We have also developed techniques for the individual resolution of ions when
the ion spacing is smaller than the beam waist. The Ramsey method technique
described in section 4.3 for the preparation of arbitrary product states works well if
the light has very little overlap on the other spins. Yet when the ion separation is of
the order of the beam waist, crosstalk significantly perturbs the other ions, reducing
95
Ins
tan
ce
of
Dis
ord
er
Ion Posi
tion
(b)
-8-6
-4-2
0
2
4
6
8
Dis
ord
er
Str
en
gth
[J max]
(a)
Figure 4.7: Fourth-order Stark shift used to generate disorder. Schematic of one par-
ticular instance of disorder (a). Each blue dot represents an ion while the black line
is a representation of the applied field. By controlling the rf frequency and ampli-
tude applied to the AOD, we can “paint” an arbitrary potential onto the ion chain.
Since these individual fields are controlled by a computer, the strength and config-
uration can be easily manipulated. Representation of all instances experimentally
realized (b). Here we show an intensity plot where each panel is a particular disorder
strength, row is a particular instance of disorder, each column is the ion position
and the color represents the strength of the potential. We experimentally realized
over 150 unique instances of disorder using the individually addressed fourth-order
Stark shift beam.
96
the fidelity of the desired state. In order to implement arbitrary state preparation
in an ion chain where the ions are closer together than the beam waist, we utilize a
different technique of single rotation known as the Rabi method (see Fig. 4.8). In
the Rabi method, we first initialize all spins in |0, 0〉 and a large fourth-order Stark
shift is applied to all ions except the ones to be rotated. We then drive a pi pulse
with a weak global carrier at the normal ωHF . All the ions with the large Stark
shift are unaffected since the carrier is off-resonant, while the unshifted ions go from
|0, 0〉 to |1, 0〉, generating the desired product state.
While both methods are useful for state preparation, only the Ramsey method
is used for analysis rotations because, in the Rabi method, the other spins undergo
a rotation around σz while the unshifted spins are flipped. This rotation could leave
an additional unwanted relative phase in the measurement, unless the Stark shift
was applied for exactly an integer 2pi time. This is difficult due to small fluctuation
in the Stark shift, making the Rabi method undesirable for analysis.
Another technique realized with this fourth-order Stark shift beam has been
the direct measurement of the Jij coupling matrix. Previously, we have recon-
structed a coupling matrix by using spectroscopic techniques [46], but with the
individual control afforded by this fourth-order Stark shift, we can now directly
measure the coupling between an arbitrary pair of ions. The procedure is similar
to the Rabi method for arbitrary state preparation. After all spins are initialized
in |0, 0〉 or |↓〉, the two spins k and l whose coupling are to be measured have a
large Stark shift applied while the rest of the spins do not. We then use a global
Raman transition to drive all the unshifted ions to one of the Zeeman levels in the
97
Ramsey evolution  
addressed ion
unaddressed ion
π pulse at unshifted
  hyperfine splitting
addressed ion
unaddressed ion
Time
AC Stark shift light
π/2 pulse π/2 pulse
final state
final state
Bloch vector evolution
Ramsey Method
Rabi Method
Figure 4.8: Schematic of individual rotation types. Upper panel is a Bloch sphere
representation of the type of rotation described in section 4.3 called the Ramsey
method. This method is ideal for both state preparation and analysis since the
phases of the unaddressed spins are unmodified. The lower panel shows what we
call the Rabi method where the unaddressed ions are shifted out of resonance while
a weak carrier rotates the unshifted spins. This method is used for state preparation
when the ion spacing is smaller than the beam waist.
98
Spin N
umberSpin Number
Sp
in-Sp
in
Co
up
lin
g(kH
z)
Short-range (α = 1.3) Long-range (α = 0.55)
Jij (kHz)
ion1
ion7
ion1
ion7 ion1
ion7
ion1
ion7
Jij (kHz)
0.0
0.5
1.0
Short-range (α = 1.3) Lo -range (α = 0.55)
Jij (kHz)
ion1
ion7
ion1
ion7 ion1
ion7
ion1
ion7
Jij (kHz)
0.0
0.5
1.0
(a) (b)
(c)
Figure 4.9: Measured Jij couplings for 7 and 10 ions. Using the individual addressing
beam, all but two spins are driven to one of the ground state manifold Zeeman states
and the two remaining spins are evolved with a M-S interaction, oscillating between
|↓↓〉 and |↑↑〉. The oscillation frequency corresponds to the interaction strength and
is used to reconstruct the coupling matrix for 7 spins with an α = 1.3 (a) and an
α = 0.55 (b). We also reconstruct a 10 ion coupling matrix with an α = 1.13 (c).
ground state manifold, namely |1,±1〉. Afterwards, a M-S interaction is applied
which induces coherent oscillations between |↓〉k |↓〉l and |↑〉k |↑〉l for only the ions
k and l since the Zeeman lines are not coupled by the M-S interaction, leaving the
rest of the ions alone. The frequency of this oscillation corresponds to the strength
of the interaction between the two ions and is extracted by fitting to the measured
populations.
This has been used to precisely measure the coupling matrix of seven ions
99
under two very different ranges, namely a short-range matrix with α = 1.3 (Fig.
4.9a) and a long-range matrix with α = 0.55 (Fig. 4.9b). The measured matrices
show good agreement with the predicted couplings and even capture the nearest-
neighbor coupling asymmetry across the chain. These Jij couplings are the very
couplings used to study a failure in quantum thermalization discussed in chapter
6. We even made a similar measurement in a system of ten ions with a range
of α ∼ 1.13 (Fig. 4.9c). This measurement highlights that the interaction graph
is really long-ranged in a ten spin system, coupling even ions 1 & 10. The long-
ranged character of the couplings is important because it prevents the Hamiltonian
discussed in [40] from being mapped to non-interacting Fermions, allowing for the
realization of a many-body localized state.
Finally, this individual addressing beam using the fourth-order Stark shift
enables dynamic individual control, allowing for the simulation of interesting systems
such as loops with non-zero magnetic flux [82], the dynamics of lattice gauge theories
[83], or even the realization of a new kind of quantum system [84]. The fourth-
order Stark shift has added a new and very powerful tool to our quantum simulator
toolbox.
100
Chapter 5: Entanglement in Large N System
As our quantum simulation system size grows, an ever present problem is the
verification and validation of the results. This problem can be broken into different
questions, one of which is, “How quantum is the system?” This question is highly
non-trivial to answer. In fact, large systems that have pursued quantum informa-
tion with a “top-down” approach are seeking to prove that their results are quantum
in origin [26, 85–87] and that these quantum processes result in a quantum advan-
tage [88, 89]. On the other hand, a “bottom-up” approach to quantum systems
begins with small systems which can be verified classically [28–30, 34–36, 90] and
seeks to increase the system size while maintaining the “quantumness” and con-
firming the dynamics. While confirming the dynamics of a system that cannot be
classically simulated is a very difficult problem, a good approximation is achieved by
measuring the imposed Hamiltonian. Verifications of an imposed Hamiltonian have
been shown that do not require an exponential number of measurements [46, 91],
ensuring that the dynamics are those believed to be applied. In conjunction with
this, a full verification of a quantum system requires that the system itself has not
undergone decoherence. One possibility to measure this decoherence is to measure
the amount of entanglement in the system, verifying that noise has not destroyed
101
this quintessential quantum phenomenon.
Entanglement of a large system is itself difficult to measure. The most di-
rect measure of entanglement is to reconstruct the density matrix of the system.
Despite density matrix reconstruction having been performed on up to 8 spins [92]
in trapped ions and small systems in photonic bits [93], complete reconstruction
quickly becomes intractable as it takes ∼ 3N measurements for an N spin system.
Thus measurement of entanglement needs to be reduced to some observable which
does not require density matrix reconstruction, but is experimentally accessible.
This has led to a variety of methods [94] proposed to measure entanglement in a
variety of systems. Some of these methods use a copy of the system to extract the
entanglement [36, 95] while others use algorithms to extract the highest probable
density matrix from an incomplete measurement set [96–98]. Witnesses constructed
for specific types of entanglement have also been successfully used in medium system
sizes [99,100], yet witnesses of this kind are generally very sensitive to experimental
noise and give only a binary answer of entangled or unknown. Ideally, in order to
act as a useful diagnostic, the measurement should somehow quantify the amount
of entanglement present in the system.
5.1 Spin Squeezing
One observable which is shown to quantify the entanglement is spin squeezing
[101,102]. Spin squeezing was originally a quantization of the metrological advantage
gained from an entangled quantum many-body system over the standard quantum
102
limit of phase resolution. Ultra-cold atomic systems have leveraged their control
in order to achieve observed metrological gains of up to 100 times the standard
quantum limit [103, 104]. This metrological gain can also be used to as a way to
quantify the k-particle entanglement [105]. Here k-particle entanglement is defined
as a state of N particles which cannot be decomposed into a convex sum of products
of density matrices with all density matrices involving less than k-particles, namely
one of the terms must be a k-particle entangled density matrix. The metrological
gain in the ultra-cold atom experiments is equivalent to the entanglement of < 1%
of the atoms in the cloud.
While quantum metrology holds great promise for high precision measurement,
spin squeezing also offers an experimentally realizable means of characterizing a
qubit system. Initial experiments in trapped ion qubits have had some success in
demonstrating metrological gain [106] as well as an entanglement of ∼ 5% of the
system. Here we will measure the spin squeezing in a large well controlled qubit
system which is used for quantum simulation.
The Hamiltonian we apply is a long-ranged transverse field Ising model,
H =
∑
i,j
Jijσ
i
xσ
j
x +B
∑
i
σiz. (5.1)
The experiment initializes all spins down, parallel to the transverse field and then
the Hamiltonian in Eq. 5.1 is applied. Since the amount of squeezing is improved
with longer range interactions, we set α = 0.74. The applied transverse field is
B = 2.5Jmax. The specific value of B is chosen to balance the sensitivity of the
103
measurement to spin depolarization due experimental noise and spin dynamics while
keeping the largest possible k-particle entanglement. By measuring the spin projec-
tions along both xˆ and zˆ, we can then construct the observables of 〈Sz〉 =
∑
i〈σiz〉
and ∆S2x =
∑
i,j〈σixσjx〉 − 〈σix〉〈σjx〉 required by the Ramsey squeezing parameter,
defined as [101]
ξ2R =
N∆S2x
〈Sz〉2 . (5.2)
In Fig. 5.1a, we show the Ramsey squeezing parameter as a function of time
for one of the largest qubit systems to date, 24 ions. The largest squeezing observed
is −2.84± 0.16 dB as compared to the intial product state. While the measurement
enhancement is not large due to the relatively small system size as compared to
5×105, it still demonstrates significant entanglement. The amount of entanglement
can be quantified using k-particle entanglement curves generated following [105].
This comparison is shown in Fig. 5.1b demonstrating 4-particle entanglement. Thus
our data is consistent with 17% of the entire system being entangled, which is much
larger than comparable experiments.
While four out of twenty-four is only a small portion of the system, numerics
only predict a maximal entanglement of 10 spins [107]. We believe that the large
disparity between the measured k-particle entanglement and the predicted value is
largely due to measurement error, particularly inter-ion crosstalk. With larger ion
numbers, the ions in the middle of the chain are pushed closer together, exacerbating
the issue of ion PSF function overlap. While this cross-talk is only a couple percent,
the variance is computed using individually resolved camera images where this ion
104
0.90 0.92 0.94 0.96 0.98 1.00
0.1
0.2
0.3
0.4
0.5
2|〈S z〉|/N
2(
Δ
S
x)
2 /
N
k-particle
Entanglement
2-particle
4-particle
6-particle
24-particle
0.0 0.1 0.2 0.3 0.4 0.5 0.60.4
0.6
0.8
1.0
1.2
Jmaxt

 R2
-2.84±0.16 dB
(a)
(b)
Figure 5.1: Measurement of spin squeezing in 24 ions. Computed Ramsey spin
squeezing parameter ξ2R from 24 spins (a). Maximum observed squeezing is -
2.84±0.16 dB. While this is small compared to ultra-cold atomic experiments, the
system size here is also significantly smaller. Also shown is ∆S2x as a function of
|〈Sz〉| (b). Curves of k-particle entanglement are generated from [105]. Data has a
maximal entanglement of four-particle entanglement. While the amount of entan-
glement is not large, theory only predicts a maximum of 10-particle entanglement
and our measurement is consistent with ∼ 17% of the entire system entangled which
is larger than comparable experiments by an order of magnitude.
105
cross-talk creates false correlations which systematically bias the variance larger
than it is. We believe that this also explains the small increase in the variance at
the second time step. This could be improved by instead using global measurements
and no longer discriminating individual spin states, but that would also sacrifice the
computational advantage of the quantum simulator. By improving our imaging
objective to a higher NA, we believe that this systematic can be eliminated even in
large chains, allowing for even deeper squeezing.
5.2 Quantum Fisher Information
Another promising witness of entanglement is the Quantum Fisher Information
(QFI) which has recently been shown to witness genuine multi-partite entanglement
[108, 109]. Genuine multi-partite entangled pure states are defined as those that
cannot be created without participation of all spins. In other words, a density
matrix ρ of N subsystems is genuinely multi-partite entangled if and only if it does
not admit any bipartite decomposition. Namely if ρ(k) is an irreducible density
matrix of k subsystems, then ρ is only genuinely multi-partite entangled if
ρ 6= ρ(k) ⊗ ρ(N−k) ∀ k. (5.3)
From a quantum metrology perspective, the QFI quantifies the sensitivity of a given
input state to a unitary transformation eiϑOˆ generated by the Hermitian operator Oˆ.
Interestingly, the QFI reveals useful entanglement even in the case of non-Gaussian
spin states where the spin squeezing witness ξR does not [110].
106
5.2.1 Directional QFI
For a pure state, the QFI for a linear local operator Oˆ reduces to the variance
of that operator:
FQ = 4(∆Oˆ)2 = 4(〈Oˆ2〉 − 〈Oˆ〉2). (5.4)
For a mixed state, the expression is more complicated, but this case will be dis-
cussed later. A direct link can now be drawn between ξR and FQ by the relation
FQ/N ≥ 1/ξ2R [110]. For states which can be fully characterized by their variance,
FQ/N ≈ 1/ξ2R but as a state becomes less Gaussian, FQ continues to indicate entan-
glement while ξR no longer does. In addition to the witness of genuine multi-partite
entanglement, there are also bounds on the QFI which determine the k-particle
entanglement present. If we write the QFI of a specific k-producible state of N
ions ρk with the observable Oˆ as FQ[ρk, Oˆ], then the QFI must satisfy the following
inequality [108,109]
FQ[ρk, Oˆ] ≤ sk2 + r2 (5.5)
where s =
⌊
N
k
⌋
is the largest integer smaller than or equal to N
k
and r = N − sk.
Thus, if the QFI violates the bound in Eq. 5.5, then there must be at least (k+ 1)-
particle entanglement. Computing this bound for k = 1, we see that FQ > N ⇒
fQ ≡ FQ/N > 1 in order for there to be at least bi-partite entanglement. It is also
clear that maximum possible value for fQ is to have N -partite entanglment in a
system of N spins, meaning that fQ ≤ N .
Let us see an example by assigning Oˆ ≡ Sz =
∑N
i σ
i
z and ρ to be a product
107
state along the xˆ direction of the Bloch sphere. If we compute the QFI, we find that
fQ[ρ, Sz] = 1, saturating the unentangled bound. If ρ was instead a product state all
along zˆ, then fQ[ρ, Sz] = 0 indicating that the direction we probe the state with the
QFI matters. This is further emphasized by instead letting ρGHZ be a GHZ state
in the zˆ direction, namely a superposition of all spins up in zˆ and all spins down
in zˆ. If we compute the QFI in each direction of the Bloch sphere, we can use the
results to generate a vector ~fQ[ρ] = {fQ[ρ, Sx], fQ[ρ, Sy], fQ[ρ, Sz]}. Applying this to
the GHZ state:
~fQ[ρGHZ ] = {1, 1, N} (5.6)
Since FQ is tied directly to metrological advantage, it is clear that there is a huge
gain in the zˆ direction while not nearly as much in the other two directions. This
implies that there is more information gained when looking in all three directions
as opposed to just one.
5.2.2 Total QFI
In fact, the bounds of k-particle entanglement can be recomputed in the con-
text of measuring in all three directions. Let a new quantity be defined which is
just the sum of the QFI in each direction of the Bloch sphere:
fQ[ρ] =
∑
l
fQ[ρ, Oˆl] (5.7)
108
where l ∈ x, y, z. The new bound for a k-producible state of N ions ρk with a linear
local observable O is [108,109]:
fQ[ρk] ≤
1
N
[
s(k2 + 2k − δk,1) + r2 + 2r − δr,1
]
(5.8)
where as before s =
⌊
N
k
⌋
is the largest integer smaller than or equal to N
k
, r = N−sk
and now δa,b is the Kronecker delta-function. For a separable state, the max value
of fQ = 2 while the maximal value of any quantum state is fQ ≤ N + 2.
Now we can revisit the example discussed above, namely a product state along
the xˆ direction where Oˆ = ~S, which gives a total QFI of fQ = 2, saturating the
separable state bound. This also gives more information about the entanglement
of the system when looking at the GHZ state in the zˆ direction ρGHZ , namely
fQ = N+2 which saturates the maximal QFI bound. While the GHZ state saturates
both fQ and fQ, there are genuinely multi-partite entangled states which are not
detected by fQ. One such example is ρDz which is a Dicke state with N/2 excitations
and is fully symmetric and an eigenstate of Jˆz. If we write out the QFI vector
~fQ[ρDz ] = {N+22 , N+22 , 0}, where none of the fQ demonstrate genuine multi-partite
entanglement, but only partial k-particle entanglement. Whereas, it is clear that if
we instead examine fQ, we find genuine multi-particle entanglement.
So far, everything stated has been in the context of pure states. When the
measured state is a mixed state, the form of the QFI itself changes. In general, for
a given density matrix ρ =
∑
k λk |k〉 〈k|, the QFI is then just FQ[ρ, Snˆ] = nˆTΓC nˆ
109
where ΓC is a matrix defined by
[ΓC ]ij = 2
∑
l,m
(λl + λm)
(
λl − λm
λl + λm
)2
〈l|Si |m〉 〈m|Sj |l〉 . (5.9)
This general statement requires knowledge of the density matrix of the state in
question, which of course takes us back to full tomography. Yet for pure states this
statement reduces to exactly [ΓC ]ij = 〈SiSj + SjSi〉/2− 〈Si〉〈Sj〉 [111] which is just
the variance of S.
Further, the QFI is convex in the state, namely fQ[p1ρ1+p2ρ2, Sl] ≤ p1fQ[ρ1, Sl]+
p2fQ[ρ2, Sl]. This implies that by computing the QFI for a pure state from a mixed
state measurement, you may be over estimating the actual size of the QFI and thus
the entanglement depth in the system. While this is true generally, there is a case
where the inequality can be rewritten as an equality. Here we define a mixed state
ρm as the following:
ρm = pρ+ (1− p) 1
2N
. (5.10)
This particular mixed state reduces to a very simple expression for the QFI by
evaluating Eq. 5.9 [109], which can be written as
fQ[ρm, Sl] =
p2
p+ (1− p)2−(N−1)fQ[ρ, Sl]. (5.11)
For large N , this expression simplifies to fQ[ρm, Sl] ∼ pfQ[ρ, Sl]. Despite this par-
ticular case, it is in general difficult to experimentally prove entanglement using
the QFI based on the variance of the observable. Some experiments have mea-
110
sured the QFI without the pure state assumption [106, 112] by using a small angle
approximation, but this technique is not generally applicable.
5.2.3 Measured QFI for 10 spins
While not proving demonstrable entanglement, the pure state QFI is still
helpful in the diagnosis of quantum simulations. Here we will consider a simulation
similar to that of [34] where the system is prepared in a product state of |↑↑ · · · ↑〉x
and then quenched to the Hamiltonian
H =
∑
i,j
Jij
2
(
σixσ
j
x + σ
i
yσ
j
y
)
. (5.12)
Previously we showed that the correlations in the xˆ direction violate the Lieb-
Robinson bounds [113] when the interaction is sufficiently long-ranged. By repeating
the quench experiment, yet also measuring in all three directions of the Bloch sphere,
we can construct the pure state QFI for the system and estimate the entanglement
in the system.
Using these measurements, we compute the total QFI fQ for a system of 10
spins (see Fig. 5.2a). The measurement shows a quick growth of entanglement in
the system, ending in a state consistent with 10-particle entanglement. Since the
QFI is convex in the state, we compared our results to exact diagonalization of Eq.
5.12.
Examination of the results showed that the dynamics were very sensitive to
small σz perturbations in the Hamiltonian. Due to small uncontrolled fourth-order
111
0 1 2 3 4 50
2
4
6
8
10
12
Jmaxt

 f Q[S
]
2 partite entanglement
6 partite entanglement
10 partite entanglement
Simple Theory
Full Theory
0
2
4
6
fQ[Sx]
0
2
4
6
fQ[Sy]
2
4
6
fQ[Sz]
(a)
(b)
Figure 5.2: Measurement of QFI in 10 Spins. Plot of 10 ion fQ as a function of time
(a). Here we see that measurements are consistent with the entirety of system being
entangled at longer times. We compare the data to theory, first with no noise (blue
line) and second with numerics which incorporate noisey Jij along with spontaneous
emission (purple line). Spontaneous emission appears to be the primary divergence
of measured dynamics from the noiseless predictions. Statistical error bars of 1 s.d.
are shown but smaller than point size. By plotting evolution of 10 ions through
the QFI vector space (b), we show that the ten ion data is consistent in all three
directions with predicted evolution (purple line) from full theory. System first starts
in product state (orange plane) and quickly evolves to violate bound for 10-particle
entanglment (green plane). Green points highlight data which violates 10-particle
bound.
112
Stark shifts arising from the Raman lasers, σz perturbations are present in the
simulation at approximately 15 Hz when the largest coupling strength is Jmax/2pi =
167 Hz. By incorporating this field into the numerics, we generate a simple model of
the system which is shown in Fig. 5.2a as the blue line. While the data is consistent
with the theory at short times, at longer times the data begins to diverge. Due to
laser intensity noise of ∼ 2%, both the strength of Jmax and the small uncontrolled
fourth-order shift are noisy. Since this noise is small and slow compared to an
experiment, we treated it perturbatively by assuming that it was constant for each
experiment. Under this assumption, we generated a Hamiltonian where Jij and the
σz were scaled by a random number chosen from a normal distribution centered
around Jmax and 15 Hz respectively. The width of the normal distribution was
set to 4%Jmax for Jij and 15 Hz for the σz. We then solved the time-dependent
Schroedinger equation for this Hamiltonian and computed the observables. This
was then repeated 100 times and the observables were then averaged together.
Finally, we incorporated the effect of spontaneous emission into the model by
assuming that the probability for an emission of a single photon per experiment
was small enough that the dynamics would be largely uneffected. This allowed us
to treat it as a measurement effect. In our system, since almost all of spontaneous
emission is Raman, e.g. emits into a different state [67], and is twice as likely to
decay into one of the Zeeman states which will project bright upon measurement,
we modified the pauli-matrices in the computation of the spin observables to be the
113
following:
σx[] =
 1− 
1 0

σy[] =
  ı˙(1− )
−ı˙ 0

σz[] =
1 
0 −(1− )

(5.13)
Using these modified Pauli-matrices, we computed the QFI where we plugged in
γ(p,|{1,0}〉) found in section 3.2.2 to the probability of a single spontaneous emission,
1− e−Nγ(p,|{1,0}〉)t.
The full theory is compared to the data in Fig. 5.2a (purple line). The data
is consistent with the full theory, showing that we have a rough characterization of
the system. We note that the most important inclusion is the spontaneous emission
probability, as this is the reason for the observed QFI decaying as a function of time.
To show that the numerics capture the observed dynamics, we also construct the
~fQ in Fig. 5.2b, where the data are the points and the theory the purple line. The
data follows the arc of the theory very well in the vector space. We highlight the
points in green which are consistent with 10-particle entanglement.
This analysis indicates that spontaneous emission is the primary cause of the
deviation of our data from the predicted values. Since spontaneous emission inde-
pendently effects all ions equally, the mixed state density matrix of our system ρobs
114
is a mixture of the coherent state density matrix ρcoh and a purely diagonal matrix
which is a mixture of all possible states ρdiag:
ρobs = pρcoh + (1− p)ρdiag (5.14)
where p is the probability of not having a spontaneous emission event at a given time.
While the exact composition of ρdiag would require reconstruction of the full density
matrix, if we approximate it as a perfectly mixed state, then the QFI reduces to Eq.
5.11. If we compare the result to the full numerics described above or the data, the
agreement is almost perfect, showing that this dissipation in our system is captured
by the QFI. While such an argument is a weak demonstration of entanglement, it
can still be used as a metric to determine the quantumness of the simulation.
5.2.4 Measured QFI for 24 spins
We emphasize this further by performing an identical experiment using 24 ions.
This system size is right on the edge of what is possible to use exact diagonalization
in order to classically simulate the system. In this case, we had our colleague
Zhexuan Gong run the numerics on a supercomputer taking a total of 2 days. By
adding only a couple more spins to the simulation, the problem becomes functionally
impossible for exact diagonalization and requires approximation methods such as
DMRG. By comparing our data to the numerics, we find that behavior is similar to
that of the 10 ion data, agreeing very well at short times but diverging at later times
(see Fig. 5.3). By applying the ideas above about the effect of spontaneous emission,
115
0 1 2 3 4 50
5
10
15
20
25
Jmaxt

 f Q[S
]
Exact Theory
Sp. Emis. Theory
0
5
10
fQ[Sx]
0
5
10
fQ[Sy]
5
10
fQ[Sz]
(a)
(b)
2 particle entanglement
11 particle entanglement
21 particle entanglement
22 particle entanglement
23 particle entanglement
24 particle entanglement
Figure 5.3: Measurement of QFI in 24 Spins. Plot of 24 ion fQ as a function of time
with pure state assumption (a). Due to size of system, exact diagonalization is not
possible on normal computer. Numerics were run on campus supercomputer, taking
2 days to complete. As few as only a couple additional spins would make such
a calculation functionally impossible, requiring approximation techniques such as
DMRG to make the problem tractable. We compare our data to the exact numerics
(blue) and also to approximation of Eq. 5.11 from spontaneous emission (purple).
We also show the 24 ion trajectory through the QFI vector space (b). While having
much of the same character as the 10 ion case, we see that the spontaneous emission
prevents the system from growing to 24-particle entanglement. Line is a guide for
the eye.
116
we compare the data to a modified theory fQobs, where fQobs = e
−Nγ(p,|{1,0}〉)t fQcoh
and fQcoh is the QFI predicted by the numerics.
The agreement between the data and the modified theory shows that the
dominant source of error in the simulation arises from spontaneous emission, which
is enhanced as we increase the system size to larger numbers. This experiment also
shows the usefulness of the QFI as a metric for how quantum a system is, despite
the fact that we do not have demonstrable entanglement. As system sizes continue
to grow, observables which are simple to measure and intimately tied to the the
coherence of a many-body system will become extremely important as a means to
verifying that the system is in some sense still quantum even when classical numerics
are no longer possible.
117
Chapter 6: Failures of Quantum Thermalization
Statistical mechanics can predict thermal equilibrium states for most classical
systems, but for an isolated quantum system there is no general understanding
on how equilibrium states dynamically emerge from the microscopic Hamiltonian
[114–120]. Exceptions to quantum thermalization have been observed, but typically
require inherent symmetries [116, 121] or noninteracting particles in the presence
of static disorder [119, 122–124]. For this reason, the predictions of many-body
localization (MBL), in which disordered quantum systems can fail to thermalize
despite strong interactions and high excitation energy, were surprising and have
attracted considerable theoretical attention [119, 120, 125–127]. This MBL phase
is predicted to emerge for a broad set of interaction ranges and disorder strengths,
though the precise phase diagram is not well known [128] since equilibrium statistical
mechanics breaks down in the MBL phase and numerical simulations are limited
to ∼ 20 particles [125, 126]. While recent experiments searching for MBL have
measured constrained mass transport and the breakdown of ergodicity in disordered
atomic systems with interactions [37,38], we have directly observed MBL in a long-
range transverse field Ising model with programmable, random disorder in a system
of ten spins [40].
118
Yet this failure of thermalization observed in MBL is only one way quantum
thermalization diverges from the classical. Despite the fact that quantum dynamics
are always unitary, we normally expect that measurements made within a subsystem
trace over the rest of the system and appear thermal because the rest of the system
acts as a thermal bath [95,129,130]. Yet, as mentioned above, there are cases where
this is not true. For integrable models an extensive number of conserved quantities
prevent the efficient exploration of phase space [116,117] and the system relaxes to
a steady-state predicted by a generalized Gibbs ensemble (GGE) [131,132] specified
by the initial values of the integrals of motion. For near-integrable systems, such
as weakly interacting ultracold gases, thermalization can still occur, but only over
extremely long time scales beyond current experimental reach [121, 133]. However,
it is possible to observe quasi-stationary states, often called prethermal, that emerge
within an experimentally accessible time scale.
Previous observations of prethermal states have focused on those described
by a GGE associated with the integrable part of the system [121, 133]. Here, we
observe a new form of prethermalization [134], where a system of interacting spins
rapidly evolves to a quasi-stationary state that cannot be described by a GGE.
This type of prethermal state arises when a system has long range interactions
and open boundaries such that the translational invariance is broken, even in the
thermodynamic limit. As a result, spin excitations feel an emergent double-well
potential whose depth grows with interaction range (Fig.6.1). Memory of the initial
state is preserved by this emergent potential, but is eventually lost due to weak
tunneling between the two wells and the interactions between spin excitations.
119
6.1 Experimental Setup
To study prethermalization in our trapped ion quantum simulator, we apply
a large transverse field Ising model Hamiltonian:
H =
∑
i<j
Jijσ
x
i σ
x
j +B
∑
i
σzi , (6.1)
where B is a uniform effective transverse field (see section 2.3.2). We set the inter-
action Jij =
Jmax
|i−j|α where Jmax/(2pi) is the maximum coupling strength which ranges
from 0.45 − 0.98 kHz, and set the magnetic field B/(2pi) = 10 kHz. We tune the
power law exponent α between α = 0.55 and 1.33 by changing the axial confinement
of the ions. With long-range interactions, H is in general non-integrable, in contrast
to the nearest-neighbor case where the 1D model is integrable [135] and thermal-
ization is anticipated in the long time limit according to eigenstate thermalization
hypothesis [136–138]. However, with a small number of spin excitations one can
map Eq. 6.1 into bosons and approximate the Hamiltonian with an integrable term
of free bosons and an integrability breaking perturbation of interacting bosons that
eventually thermalizes the system at times longer than ∼ B2/J3 (see section 6.3).
We initialize the chain of seven spins by optically pumping all spins to the
|↓z〉 state, and then use a tightly focused individual-ion addressing laser to excite
a single spin on one end of the chain to the |↑z〉 state as seen in Fig. 6.1a (see
chapter 4). The spins then evolve under Eq. 6.1 and we measure the time evolu-
tion of the spin projection in the zˆ-basis. For the shortest range interactions we
120
Small α
Large α
Short range
(b)(a)
(d)(c)
Ueff
Ueff
Initial spin
exitation Generalized Gibbs
 Ensemble (GGE) prediction
GGE prediction
Long range
ΔE/Jmax
Ueff
α
1     2     3     4     5     6     7Ion # 1     2     3     4     5     6     7Ion #
1     2     3     4     5     6     7Ion #
E
xc
ita
tio
n 
pr
ob
ab
ili
ty
E
xc
ita
tio
n 
pr
ob
ab
ili
ty
E
xc
ita
tio
n 
pr
ob
ab
ili
tym|ψ0         ψ0|n         | |
2
α=1.33
α=0.55
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.1
0.2
0.3
0.4
0.5
0.0
0.1
0.2
0.3
0.4
0.5
Figure 6.1: Emergent double well potential. (a) The spin chain starts with a single
spin excitation on the left end. (b) For short range interactions (α = 1.33), we map
the system to a particle in a 1D square well where the excitation becomes symmet-
rically distributed across the chain as predicted by the GGE, 〈σzi 〉GGE. (d) However,
for long-range interactions (α = 0.55) the Hamiltonian resembles a double-well po-
tential which prevents the efficient transfer of the spin and the excitation location
retains memory of the initial state and resides primarily on the left half, in con-
trast to 〈σzi 〉GGE. (c) (Double well gives rise to ) near-degenerate eigenstates which
emerge as α is decreased as seen in the calculated energy difference between all pairs
of eigenstates versus α with amplitude weighted by the product of the eigenstates
overlap with the initial state . Because of this near degeneracy, the initial phase
relationship between these two eigenstates remains frozen for the experimentally
accessible timescales.
121
realize (α = 1.33) the system rapidly evolves to a prethermal state predicted by the
GGE associated with the integrals of motion corresponding to the momentum space
occupation number of the non-interacting bosons (see section 6.3) (Fig. 6.1b).
However, in the long-range interacting case (α = 0.55) we see the location of
the spin excitation reaches an equilibrium value that retains a memory of the initial
state (Fig. 6.1d) out to the longest experimentally achievable time of 25/Jmax. This
prethermal state is in obvious disagreement with both a thermal state and the GGE
prediction, as both maintain the right-left symmetry of the system.
In the thermodynamic limit, the spin-wave boson model for short-range inter-
actions are similar to those of a free-particle in a square-well. However, for long-range
interactions, there is an emergent double-well potential and the spin-wave model has
an extensive number of near-degenerate eigenstates. These eigenstates, which are
superpositions of its symmetric and antisymmetric ground states, correspond to spin
excitations in the left and right potential well. For seven spins, we calculate the en-
ergy difference between all pairs of eigenstates for seven spins and plot them with
respect to α in Fig. 6.1c where the amplitudes, |〈m|ψ0〉〈ψ0|n〉|2, are the product of
each eigenstates’, |m〉, |n〉, overlap with the initial state. With α = 0.55 the two
lowest energy states are almost degenerate, with an energy difference one-thousand
times smaller than Jmax, due to the tunneling rate which is exponentially small in
the height of the double well. As a result the spin excitation will remain in its initial
well until it tunnels into the other well at a much longer time.
122
G
G
E
Jmax t
G
G
E
G
G
E
Jmax t
G
G
E Spin index
Spin 
index
Spin 
index
Spin 
index
Jmax t Jmax t
Jmax tJmax t Jmax t
Short-range Long-range
S
in
gl
e 
sp
in
-fl
ip ψ0         =
ψ0         =|
|
ψ0         =
ψ0         =|
| ψ0         =
ψ0         =|
|
D
ou
bl
e 
sp
in
-fl
ip
ψ0         =
ψ0         =|
|
Figure 6.2: Measured location of spin excitation. The average position of the spin
excitation, 〈C〉, is plotted for an initial excitation on the left (dark blue) and right
(dark red) along with their cumulative time average (blue and red), for short-range
(α = 1.33) and long-range (α = 0.55) interactions with initial states with one
spin excitation (left panel) and two spin excitations (right panel). The cumulative
time average of the deviation of the postselected individual spin projections from
the generalized Gibbs ensemble, 〈σzi 〉 − 〈σzi 〉GGE, is also plotted. For short-range
interactions the spins quickly thermalize to a value predicted by the GGE, but for
the long-range interacting system a prethermal state emerges.
6.2 Observation of Prethermalization
To better characterize the “location” of the spin excitation we construct the
observable
C =
∑
i
2i−N − 1
N − 1
σzi + 1
2
, (6.2)
where N is the number of ions. The expectation value of C varies between -1 and
1 for a spin excitation on the left and right ends respectively. Due to the spatial
inversion symmetry of the spin-wave boson model and Eq. 6.1, the GGE predicted
and thermal values of 〈C〉 are zero. In Fig. 6.2, we plot 〈C〉 along with its cumulative
time average 〈C〉 for α = 1.33 (short-range) and α = 0.55 (long-range). To further
emphasize the asymmetry of the prethermalization we prepare initial states with a
123
single-spin excitation on the left or right ends of the chain.
With short-range interactions the system quickly relaxes to the prethermal
value predicted by the GGE where 〈C〉 is near zero irrespective of the initial state.
〈Ct=50/Jmax〉 = −0.018(4) (−0.032(4)) for the initial spin excitation on the right
(left). But with long-range interactions, the prethermal state that emerges retains
a clear memory of the initial conditions as 〈Ct=25/Jmax〉 = 0.205(4) (-0.179(4)) for
the excitation starting on the right (left), showing a clear discrepancy with both the
thermal and GGE predictions.
It is useful to talk about thermalization in terms of local quantities since
subsystems must use the rest of the chain as a heat bath. In Fig. 6.2 we plot the
cumulative time average of the deviation of the single spin magnetizations 〈σzi 〉 from
the GGE. Here we postselect the data for the correct number of spin excitations to
eliminate the effects of imperfect state preparation and small deviations from our
model Hamiltonian due to unwanted excitations of the phonon modes. For the
short-range interactions we see the cumulative time average quickly converge to
the GGE predicted value. However, with the long-range interactions we see the
individual spins reach a steady state that does not match the GGE as the pair of
near degenerate eigenstates prevents the efficient transfer of the excitation across
the chain.
We show the prethermal state’s robustness to weak interactions, similar to
MBL [119,120], by preparing initial states with two spin excitations. In this case, the
multiple spin flips increase the size of the integrability breaking part of the Hamil-
tonian which represents weak interactions between the spin-wave bosons. How-
124
Jmax t
(a)
(b)
Figure 6.3: Comparison of prethermal data to exact numerics. Comparison between
a numerical simulation of the transverse field Ising model and experimental data
for a single spin excitation (top panel) and a double excitation (bottom panel).
The inset shows excellent agreement between experiment and numerics accounting
for experimental noise for 〈C〉. For both the single and double spin flip cases,
the prethermal state persists for much longer than the experimental timescale and
eventually relaxes to the thermal state.
125
ever, the prethermal state still persists. For better contrast between the prethermal
and GGE predicted values of 〈C〉, we flip the second and fourth spin such that
|ψ0〉 = |↓↑↓↑↓↓↓〉. We emphasize that the observed prethermalization and depar-
ture from the GGE prediction is not sensitive to the specific choice of initial state.
But for a small sized system, these initial states offer us the maximum signal to noise
ratio. As before, we also prepare the mirror image of the initial state by exciting
the fourth and sixth ions (ψ0〉 = |↓↓↓↑↓↑↓〉). We observe relaxation to the value
predicted by the GGE for short-range interactions, but with long-range interactions
we see a prethermal state that strongly deviates from the GGE (bottom panel of
Fig. 6.2).
In Fig. 6.3 we plot the experimental evolution of the prethermal state for both
the double and single spin flip initial states along with a long-time numerical simu-
lation under Eq. 6.1 with an inset comparing the short-time dynamics to a model
accounting for known experimental noise. There is excellent agreement between the
two and confirmation that the prethermal states persist well beyond the current
experimental time limit. Due to the interactions between the spin-wave bosons, the
system will eventually relax to the thermal equilibrium in the thermodynamic limit,
however relaxation to the GGE may or may not be seen depending on the range of
interactions.
To demonstrate that the prethermal state we observe is not due to small system
size, we repeat the experiment in a chain of 22 ions, the largest ion chain used for
quantum simulation in the literature to date. The time evolution and cumulative
time average of 〈C〉 are depicted in Fig. 6.4. Experimentally, as well as in the
126
0 6 12 18 24 30 36
-1.0
-0.5
0.0
Jmax t
C
0.5
1.0
22 spins all up
Single spin excitation
End of spin evolution
Figure 6.4: Scaling to large system size. Time evolution (light blue) and cumulative
time average (orange) of 〈C〉 with false color pictures of the 22 ion chain, where the
brightness of each ion is determined by the value of 〈σzi 〉. The ions fluoresce during
detection when in the |↑〉 state (top picture). We initialize the chain with a single
spin excitation on the left end of the chain (middle picture). After evolving for 36
Jmax the spin excitation is delocalized, but its average position remains on the left
half of the chain (bottom picture).
analytic result in the thermodynamic limit, we see the system relaxes to a similar
quasi-equilibrium state as before that clearly has memory of the initial state.
We point out the observed prethermalization and deviation from the GGE
should disappear if the system is subject to periodic boundary conditions, regard-
less of system size. Here we emphasize that the long-range interactions make the
boundary conditions relevant for bulk properties and thus changing the boundary
conditions can impact an extensive number of eigenstates. The effects of long-
range interactions on many-body dynamics are, however, far richer than the effect
discussed in the current experiment. For sufficiently long-range interactions, the no-
tion of locality breaks down and quasi-particles in the system can travel at divergent
127
velocities for thermodynamic systems, potentially leading to dramatically different
thermalization/prethermalization time scales in certain systems [34,35,139].
6.3 Spin-boson mapping and the GGE
To explain our observed prethermalization, it is convenient to map the spins
into bosons by using the Holstein-Primakoff transformation: σzi = 2a
†
iai − 1, σ+i =
a†i
√
1− a†iai. We will assume that the average spin excitation density n¯ =
∑
i〈a†iai〉/N
is much smaller than 1. This assumption is justified in our experiments because (1)
our initial states have small spin excitation densities and (2) we set max(Jij)  B
so the amount of n¯ that will be dynamically created is small [∼ (max(Jij)/B)2].
Therefore to the lowest order we can approximate σ+i ≈ a†i , and Eq. 6.1 reduces to
an integrable Hamiltonian H0 made of non-interacting bosons
H0 =
∑
i<j
Jij(a
†
iaj + a
†
ia
†
j + h.c.) + 2B
∑
i
a†iai
H1 = H −H0.
(6.3)
Here H1 contains interactions between the bosons which are parametrically small in
n¯, and, as a result, we can treat H1 as a perturbation to H0. Thus, it is natural to
expect the system (in the thermodynamic limit) to first relax to a prethermal state
described by the GGE of H0, and to later relax to a thermal state described by
the full H. Naively, we expect the thermalization to happen at a time scale much
longer than the relaxation to GGE, based on the different energy scales of H0 and
H1. However, this is not always the case, as discussed below.
128
Integrable models have many conserved quantities that are not taken into
account by traditional ensembles from statistical mechanics. Thus, the GGE was
developed to make predictions about equilibrium values of observables in these sys-
tems by incorporating the additional integrals of motion [131, 132, 140–143]. The
GGE predicted density matrix is given by:
ρGGE =
∑
k
e−λk Iˆk/tr(e−λk Iˆk) (6.4)
where λk is the overlap of the initial state with the integral of motion, Iˆk. We observe
relaxation to the GGE if for any local observable Oˆ, 〈O(t)〉 ≈ tr(OρGGE). The GGE
is equivalent to the grand canonical ensemble for a general quantum system when
the only conserved quantities are particle number and total energy [131].
To explicitly define the GGE of H0, we would need to first diagonalize H0 and
find the integrals of motion. Note that Eq. 6.1 only involves Jij for i < j. For
convenience, we will define a new matrix J such that Jii = 0 and Jij = Jji = Jij
for i < j. H0 can be rewritten as
H0 =
∑
i,j
Jij
[
a†iaj +
1
2
(a†ia
†
j + aiaj)
]
+ 2B
∑
i
a†iai. (6.5)
An orthogonal matrix V can be used to diagonalize the matrix J as∑i,j VikJijVjk′ =
νkδkk′ , where {νk} are the eigenvalues of matrix J . Introducing ck =
∑
i Vikai, we
have
H0 =
∑
k
[
(νk + 2B)c
†
kck +
1
2
νk(c
†
kc
†
k + ckck)
]
. (6.6)
129
Next, we perform a Bogoliubov transformation ck = cosh(θk)dk − sinh(θk)d†k with
θk =
1
2
tanh−1( νk
νk+2B
) to fully diagonalize H0:
H0 =
∑
k
kd
†
kdk k ≡ 2
√
B(B + νk), (6.7)
The GGE for H0 is defined as
ρGGE =
e−
∑
k λkd
†
kdk
Tr(e−
∑
k λkd
†
kdk)
, (6.8)
with λk determined by 〈d†kdk〉0 = 〈d†kdk〉GGE, where the notation 〈· · · 〉0 denotes the
expectation value in the initial state |ψ0〉, and the notation 〈· · · 〉GGE denotes the
expectation value in ρGGE. Using the fact that our initial state |ψ0〉 is always a Fock
state in the basis of {a†iai}, the value of 〈d†kdk〉0 can be calculated by the following
formula
〈d†kdk〉0 = cosh(2θk)
∑
i
V2ik〈a†iai〉0 + sinh2(θk). (6.9)
To calculate the expectation values of σzi = 2a
†
iai − 1 in GGE, we use the
following equations
〈a†iai〉GGE = 〈
∑
k,k′
VikVik′c†kck′〉GGE,
=
∑
k
V2ik〈cosh(2θk)d†kdk + sinh2(θk)〉GGE,
=
∑
k
V2ik[cosh(2θk)〈d†kdk〉0 + sinh2(θk)],
(6.10)
130
where we use the fact that ρGGE is diagonal in the Fock basis of {d†kdk}, so 〈d†kdk′〉GGE =
〈d†kd†k′〉GGE = 0 for k 6= k′.
6.4 Origin of the double well
The low-energy eigenstates of a massive particle in a double-well-shaped po-
tential are localized orbitals inside either well (Fig. 6.1b). The tunneling rate (which
is exponentially small in the height of the barrier) splits the degeneracy of the lo-
calized orbitals in each well and leads to pairs of symmetric and antisymmetric
wavefunctions (upon the spatial inversion of the chain). This physical picture ex-
plains the observed prethermalization: Initial excitations placed in the left half of
the chain will be localized for an extended period of time (inversely proportional to
the tunneling rate) under the evolution of H0, until the tunneling between the two
wells eventually delocalizes the excitations.
Interestingly, the double-well-shaped potential is emergent in our system, be-
cause the non-interacting Hamiltonian H0 contains no inhomogeneous potential.
This emergent inhomogeneity is somewhat a surprising effect because the motional
Hamiltonian of the ions in the trap and the spin-motion couplings induced by the
Raman lasers are all homogeneous. The key is that the long-range interactions make
our spin Hamiltonian Eq. 6.1 no longer translationally invariant (even in the ther-
modynamic limit). This is in contrast to an open-boundary spin chain with only
short-range interactions, where it is usually safe to assume translational invariance
for sufficiently large system sizes.
131
The notion of boundary starts to break down for sufficiently long-ranged in-
teractions, and therefore we cannot attribute the observed prethermalization to
boundary effects. Usually, boundary effects only affect a finite number of eigen-
states and do not affect local quenches in the bulk. However, here there is extensive
number of eigenstates that are localized in the two wells described by the potential
U(z), and excitations placed with an extensive number of lattice sites away from
the edges are still subject to qualitatively similar dynamics.
Finally, we point out that interactions in H1 can also delocalize the initial
spin excitations placed in one of the wells, and eventually thermalize the system.
As a result, there is an interesting interplay between the timescale of the prether-
malization to GGE and the timescale of thermalization. If the interactions in H1
are sufficiently weaker than the kinetic tunneling rate in H0 (which can be achieved
by turning down the magnetic field strength B, or changing the range of interac-
tions), then we expect the system to have two prethermal phases, with the observed
prethermalization followed by prethermalization to the GGE of H0, and further
followed by the thermalization process. On the contrary, if the tunneling rate is
sufficiently smaller than the interaction strength of H1, then the prethermal states
described by the GGE of H0 may never appear during the time evolution, and we
will find thermalization directly following the observed prethermalization. These
interesting multi-stage relaxation processes will require future experimental investi-
gations, where longer coherence time and larger number of the spins are technically
achievable.
132
Bibliography
[1] P. Shor, “Algorithms for quantum computation: Discrete logarithms and fac-
toring,” Proceedings of the Annual Symposium on Foundations of Computer
Science, pp. 124–134, 1994.
[2] L. K. Grover, “Quantum mechanics helps in searching for a needle in a
haystack,” Phys. Rev. Lett., vol. 79, p. 325, 1997.
[3] M. Mosca, Quantum Algorithms, pp. 7088–7118. New York, NY: Springer
New York, 2009.
[4] A. M. Childs and W. van Dam, “Quantum algorithms for algebraic problems,”
Rev. Mod. Phys., vol. 82, pp. 1–52, Jan 2010.
[5] B. M. Terhal, “Quantum error correction for quantum memories,” Rev. Mod.
Phys., vol. 87, pp. 307–346, Apr 2015.
[6] T. D. Ladd, F. Jelezko, R. Laflamme, Y. Nakamura, C. Monroe, and J. L.
O’Brien, “Quantum computers,” Nature, vol. 464, p. 45, 2010.
[7] I. Bloch, J. Dalibard, and W. Zwerger, “Many-body physics with ultracold
gases,” Rev. Mod. Phys., vol. 80, p. 885, 2008.
[8] F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. L. Hollen-
berg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, “Silicon
quantum electronics,” Rev. Mod. Phys., vol. 85, pp. 961–1019, Jul 2013.
[9] J. Clarke and F. K. Wilhelm, “Superconducting quantum bits,” Nature,
vol. 453, pp. 1031–1042, June 2008.
[10] R. Blatt and D. Wineland, “Entangled states of trapped atomic ions,” Nature,
vol. 453, p. 1008, 2008.
[11] C. J. Ballance, T. P. Harty, N. M. Linke, M. A. Sepiol, and D. M. Lucas,
“Laser-driven quantum logic gates with precision beyond the fault-tolerant
threshold.” arXiv:1512.04600, 2015.
133
[12] J. P. Gaebler, T. R. Tan, Y. W. Y. Lin, R. Bowler, A. C. Keith, S. Glancy,
K. Coakley, E. Knill, D. Leibfried, and D. J. Wineland, “High-fidelity universal
gate set for 9be+ ion qubits.” arXiv:1604.00032, 2016.
[13] Z. Chen, J. Kelly, C. Quintana, R. Barends, B. Campbell, Y. Chen, B. Chiaro,
A. Dunsworth, A. G. Fowler, E. Lucero, E. Jeffrey, A. Megrant, J. Mutus,
M. Neeley, C. Neill, P. J. J. O’Malley, P. Roushan, D. Sank, A. Vainsencher,
J. Wenner, T. C. White, A. N. Korotkov, and J. M. Martinis, “Measuring and
suppressing quantum state leakage in a superconducting qubit,” Phys. Rev.
Lett., vol. 116, p. 020501, Jan 2016.
[14] R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C.
White, J. Mutus, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro,
A. Dunsworth, C. Neill, P. O/’Malley, P. Roushan, A. Vainsencher, J. Wen-
ner, A. N. Korotkov, A. N. Cleland, and J. M. Martinis, “Superconducting
quantum circuits at the surface code threshold for fault tolerance,” Nature,
vol. 508, pp. 500–503, Apr. 2014.
[15] T. Monz, D. Nigg, E. A. Martinez, M. F. Brandl, P. Schindler, R. Rines, S. X.
Wang, I. L. Chuang, and R. Blatt, “Realization of a scalable shor algorithm,”
Science, vol. 351, no. 6277, pp. 1068–1070, 2016.
[16] S. Debnath, N. M. Linke, C. Figgatt, K. A. Landsman, K. Wright, and
C. Monroe, “Demonstration of a programmable quantum computer module.”
arXiv:1603.04512, March 2016.
[17] S. Lloyd, “Universal quantum simulators,” Science, vol. 273, p. 1073, 1996.
[18] R. Feynman, “Simulating physics with computers,” International Journal of
Theoretical Physics, vol. 21, p. 467, 1982.
[19] I. M. Georgescu, S. Ashhab, and F. Nori, “Quantum simulation,” Rev. Mod.
Phys., vol. 86, pp. 153–185, Mar 2014.
[20] I. Bloch, J. Dalibard, and S. Nascimbene, “Quantum simulations with ultra-
cold gases,” Nature Physics, vol. 8, p. 267, 2012.
[21] R. Blatt and C. F. Roos, “Quantum simulations with trapped ions,” Nature
Physics, vol. 8, p. 277, 2012.
[22] A. Aspuru-Guzik and P. Walther, “Photonic quantum simulators,” Nature
Physics, vol. 8, p. 285, 2012.
[23] A. A. Houck, H. E. Tureci, and J. Koch, “On-chip quantum simulation with
superconducting circuits,” Nat Phys, vol. 8, pp. 292–299, Apr. 2012.
[24] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda,
“A quantum adiabatic evolution algorithm applied to random instances of an
np-complete problem,” Science, vol. 292, p. 472, 2001.
134
[25] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson,
R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud,
J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov,
C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang,
B. Wilson, and G. Rose, “Quantum annealing with manufactured spins,” Na-
ture, vol. 473, pp. 194–198, May 2011.
[26] S. Boixo, T. F. Ronnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M.
Martinis, and M. Troyer, “Evidence for quantum annealing with more than
one hundred qubits,” Nat Phys, vol. 10, pp. 218–224, Mar. 2014.
[27] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, and M. Greiner,
“Quantum simulation of antiferromagnetic spin chains in an optical lattice,”
Nature, vol. 472, p. 307, 2011.
[28] R. Islam, C. Senko, W. C. Campbell, S. Korenblit, J. Smith, A. Lee, E. E.
Edwards, C.-C. J. Wang, J. K. Freericks, and C. Monroe, “Emergence and
frustration of magnetic order with variable-range interactions in a trapped ion
quantum simulator,” Science, vol. 340, p. 583, 2013.
[29] P. Richerme, C. Senko, J. Smith, A. Lee, S. Korenblit, and C. Monroe, “Exper-
imental performance of a quantum simulator: Optimizing adiabatic evolution
and identifying many-body ground states,” Phys. Rev. A, vol. 88, p. 012334,
2013.
[30] R. Barends, A. Shabani, L. Lamata, J. Kelly, A. Mezzacapo, U. L. Heras,
R. Babbush, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro,
A. Dunsworth, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley,
C. Neill, P. J. J. OMalley, C. Quintana, P. Roushan, D. Sank, A. Vainsencher,
J. Wenner, T. C. White, E. Solano, H. Neven, and J. M. Martinis, “Digi-
tized adiabatic quantum computing with a superconducting circuit,” Nature,
vol. 534, pp. 222–226, June 2016.
[31] T. Uehlinger, G. Jotzu, M. Messer, D. Greif, W. Hofstetter, U. Bissbort, and
T. Esslinger, “Artificial graphene with tunable interactions,” Phys. Rev. Lett.,
vol. 111, p. 185307, Oct 2013.
[32] N. Goldman, G. Juzelinas, P. hberg, and I. B. Spielman, “Light-induced gauge
fields for ultracold atoms,” Reports on Progress in Physics, vol. 77, no. 12,
p. 126401, 2014.
[33] M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P. Schauss, T. Fukuhara,
C. Gross, I. Bloch, C. Kollath, and S. Kuhr, “Light-cone-like spreading of
correlations in a quantum many-body system,” Nature, vol. 481, p. 484, 2012.
[34] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Micha-
lakis, A. V. Gorshkov, and C. Monroe, “Non-local propagation of correlations
135
in quantum systems with long-range interactions,” Nature, vol. 511, p. 198,
2014.
[35] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt, and C. F.
Roos, “Quasiparticle engineering and entanglement propagation in a quantum
many-body system,” Nature, vol. 511, p. 202, 2014.
[36] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and
M. Greiner, “Measuring entanglement entropy in a quantum many-body sys-
tem,” Nature, vol. 528, pp. 77–83, Dec. 2015.
[37] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lschen, M. H. Fischer, R. Vosk,
E. Altman, U. Schneider, and I. Bloch, “Observation of many-body local-
ization of interacting fermions in a quasi-random optical lattice,” Science,
vol. 349, p. 842, 2015.
[38] S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco, “Disorder-induced
localization in a strongly correlated atomic hubbard gas,” Phys. Rev. Lett.,
vol. 114, p. 083002, Feb 2015.
[39] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khe-
mani, D. A. Huse, I. Bloch, and C. Gross, “Exploring the many-body localiza-
tion transition in two dimensions,” Science, vol. 352, no. 6293, pp. 1547–1552,
2016.
[40] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl,
D. A. Huse, and C. Monroe, “Many-body localization in a quantum simu-
lator with programmable random disorder,” Nat Phys, vol. advance online
publication, pp. –, June 2016.
[41] J. P. Covey, S. A. Moses, M. Garttner, A. Safavi-Naini, M. T. Miecnikowski,
Z. Fu, J. Schachenmayer, P. S. Julienne, A. M. Rey, D. S. Jin, and J. Ye,
“Doublon dynamics and polar molecule production in an optical lattice,” Nat
Commun, vol. 7, pp. –, Apr. 2016.
[42] M. Saffman, T. G. Walker, and K. Mølmer, “Quantum information with ryd-
berg atoms,” Rev. Mod. Phys., vol. 82, pp. 2313–2363, Aug 2010.
[43] R. Loew, H. Weimer, J. Nipper, J. B. Balewski, B. Butscher, H. P. Buechler,
and T. Pfau, “An experimental and theoretical guide to strongly interacting
rydberg gases,” Journal of Physics B: Atomic, Molecular and Optical Physics,
vol. 45, no. 11, p. 113001, 2012.
[44] O. Firstenberg, T. Peyronel, Q.-Y. Liang, A. V. Gorshkov, M. D. Lukin, and
V. Vuletic, “Attractive photons in a quantum nonlinear medium,” Nature,
vol. 502, p. 71, 2013.
136
[45] K. Kim, S. Korenblit, R. Islam, E. E. Edwards, M.-S. Chang, C. Noh,
H. Carmichael, G.-D. Lin, L.-M. Duan, C.-C. J. Wang, J. K. Freericks, and
C. Monroe, “Quantum simulation of the transverse field ising model with
trapped ions,” New J. Physics, vol. 13, p. 105003, 2011.
[46] C. Senko, J. Smith, P. Richerme, A. Lee, W. C. Campbell, and C. Monroe,
“Coherent imaging spectroscopy of a quantum many-body spin system,” Sci-
ence, vol. 345, p. 430, 2014.
[47] B. Neyenhuis, J. Smith, A. C. Lee, J. Zhang, P. Richerme, P. W. Hess, Z.-X.
Gong, A. V. Gorshkov, and C. Monroe, “Observation of prethermalization in
long-range interacting spin chains.” In preparation.
[48] K. R. Islam, Quantum Simulation of Interacting Spin Models with Trapped
Ions. PhD thesis, University of Maryland, College Park, 2012.
[49] S. Korenblit, Quantum Simulations of the Ising Model: Devil’s Staircase and
Arbitrary Lattice Proposal. PhD thesis, University of Maryland, College Park,
2013.
[50] C. Senko, Dynamics and excited states of quantum many-body spin chains with
trappedions. PhD thesis, University of Maryland, 2014.
[51] D. Hayes, D. N. Matsukevich, P. Maunz, D. Hucul, Q. Quraishi, S. Olmschenk,
W. Campbell, J. Mizrahi, C. Senko, and C. Monroe, “Entanglement of atomic
qubits using an optical frequency comb,” Phys. Rev. Lett., vol. 104, p. 140501,
2010.
[52] D. J. Wineland, C. Monroe, W. M. Itano, D. Leibfried, B. E. King, and
D. M. Meekhof, “Experimental issues in coherent quantum-state manipula-
tion of trapped atomic ions,” Journal of Research of the National Institute of
Standards and Technology, vol. 103, p. 259, 1998.
[53] E. T. Jaynes and F. W. Cummings, “Comparison of quantum and semiclassical
radiation theories with application to the beam maser,” Proc. IEEE, vol. 51,
p. 89, 1963.
[54] K. Molmer and A. Sorensen, “Multiparticle entanglement of hot trapped ions,”
Phys. Rev. Lett., vol. 82, p. 1835, 1999.
[55] E. Solano, R. L. de Matos Filho, and N. Zagury, “Deterministic bell states
and measurement of the motional state of two trapped ions,” Phys. Rev. A,
vol. 59, p. R2539, 1999.
[56] S.-L. Zhu, C. Monroe, and L.-M. Duan, “Trapped ion quantum computation
with transverse phonon modes,” Phys. Rev. Lett., vol. 97, p. 050505, 2006.
137
[57] T. Choi, S. Debnath, T. A. Manning, C. Figgatt, Z.-X. Gong, L.-M. Duan,
and C. Monroe, “Optimal quantum control of multimode couplings between
trapped ion qubits for scalable entanglement,” Phys. Rev. Lett., vol. 112,
p. 190502, May 2014.
[58] P. J. Lee, K.-A. Brickman, L. Deslauriers, P. C. Haljan, L.-M. Duan, and
C. Monroe, “Phase control of trapped ion quantum gates,” Journal of Optics
B, vol. 7, p. S371, 2005.
[59] H. F. Trotter, “On the product of semi-groups of operators,” Proceedings of
the American Mathematical Society, vol. 10, p. 545, 1959.
[60] B. P. Lanyon, C. Hempel, D. Nigg, M. Mu¨ller, R. Gerritsma, F. Za¨hringer,
P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, M. Hennrich,
P. Zoller, R. Blatt, and C. F. Roos, “Universal digital quantum simulation
with trapped ions,” Science, vol. 334, no. 6052, pp. 57–61, 2011.
[61] C. Balzer, A. Braun, T. Hannemann, C. Paape, M. Ettler, W. Neuhauser, and
C. Wunderlich, “Electrodynamically trapped yb+ ions for quantum informa-
tion processing,” Phys. Rev. A, vol. 73, p. 041407, Apr 2006.
[62] S. Olmschenk, K. C. Younge, D. L. Moehring, D. N. Matsukevich, P. Maunz,
and C. Monroe, “Manipulation and detection of a trapped yb+ hyperfine
qubit,” Phys. Rev. A, vol. 76, p. 052314, 2007.
[63] R. Islam, W. C. Campbell, T. Choi, S. M. Clark, S. Debnath, E. E. Edwards,
B. Fields, D. Hayes, D. Hucul, I. V. Inlek, K. G. Johnson, S. Korenblit, A. Lee,
K. W. Lee, T. A. Manning, D. N. Matsukevich, J. Mizrahi, Q. Quraishi,
C. Senko, J. Smith, and C. Monroe, “Beat note stabilization of mode-locked
lasers for quantum information processing,” Optics Letters, vol. 39, p. 3238,
2014.
[64] W. C. Campbell, J. Mizrahi, Q. Quraishi, C. Senko, D. Hayes, D. Hucul,
D. N. Matsukevich, P. Maunz, and C. Monroe, “Ultrafast gates for single
atomic qubits,” Phys. Rev. Lett., vol. 105, p. 090502, 2010.
[65] J. Mizrahi, Ultrafast control of spin and motion in trapped Ions. PhD thesis,
University of Maryland, 2013.
[66] D. J. Wineland, M. Barrett, J. Britton, J. Chiaverini, B. DeMarco, W. M.
Itano, B. Jelenkovic, C. Langer, D. Leibfried, V. Meyer, T. Rosenband, and
T. Schaetz, “Quantum information processing with trapped ions,” Philosoph-
ical Transactions of the Royal Society A, vol. 361, pp. 1349–1361, 2003.
[67] R. Ozeri, W. M. Itano, R. B. Blakestad, J. Britton, J. Chiaverini, J. D. Jost,
C. Langer, D. Leibfried, R. Reichle, S. Seidelin, J. H. Wesenberg, and D. J.
Wineland, “Errors in trapped-ion quantum gates due to spontaneous photon
scattering,” Phys. Rev. A, vol. 75, p. 042329, 2007.
138
[68] H. J. Metcalf and P. van der Straten, Laser Cooling and Trapping. Springer-
Verlag, 1999.
[69] A. Kramida, Yu. Ralchenko, J. Reader, and and NIST ASD Team.
NIST Atomic Spectra Database (ver. 5.3), [Online]. Available:
http://physics.nist.gov/asd [2016, March 31]. National Institute of
Standards and Technology, Gaithersburg, MD., 2015.
[70] B. C. Fawcett and M. Wilson, “Computed oscillator strengths, Lande´ g val-
ues, and lifetimes in Yb II,” Atomic Data and Nuclear Data Tables, vol. 47,
pp. 241–317, 1991.
[71] E. Bie´mont, J.-F. Dutrieux, I. Martin, and P. Quinet, “Lifetime calculations
in Yb II,” J. Phys. B: At. Mol. Opt. Phys., vol. 31, pp. 3321–3333, 1998.
[72] J. Hynecek and T. Nishiwaki, “Excess noise and other important characteris-
tics of low light level imaging using charge multiplying ccds,” IEEE Transac-
tions on Electron Devices, vol. 50, pp. 239–245, Jan 2003.
[73] A. C. Lee, “Ytterbium ion qubit state detection on an iccd camera,” b.sc.
thesis, University of Maryland, College Park, 2012.
[74] C. Shen and L.-M. Duan, “Correcting detection errors in quantum state engi-
neering through data processing,” New Journal of Physics, vol. 14, p. 053053,
2012.
[75] d. Lichtman, Martin Tom, Coherent operations, entanglement, and progress
toward quantum search in a large 2D array of neutral atom qubits. [Madison,
Wis.] : [University of Wisconsin–Madison], 2015., 2015.
[76] J. J. Bollinger, D. J. Heinzen, W. M. Itano, S. . L. Gilbert, and D. J. Wineland,
“A 303-mhz frequency standard based on trapped be ions,” IEEE Transactions
on Instrumentation and Measurement, vol. 40, no. 2, pp. 126–128, 1991.
[77] P. Fisk, M. Sellars, M. Lawn, and C. Coles, “Accurate measurement of the
12.6 ghz ”clock” transition in trapped /sup 171/yb/sup +/ ions,” Ultrasonics,
Ferroelectrics, and Frequency Control, IEEE Transactions on, vol. 44, pp. 344–
354, March 1997.
[78] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Infor-
mation: 10th Anniversary Edition. New York, NY, USA: Cambridge Univer-
sity Press, 10th ed., 2011.
[79] J. Mizrahi, B. Neyenhuis, K. G. Johnson, W. C. Campbell, C. Senko, D. Hayes,
and C. Monroe, “Quantum control of qubits and atomic motion using ultrafast
laser pulses,” Applied Physics B, vol. 114, no. 1, pp. 45–61, 2013.
[80] Zemax LLC, “Zemax (version 13.2),” 2013.
139
[81] N. F. Ramsey, “Experiments with separated oscillatory fields and hydrogen
masers,” Rev. Mod. Phys., vol. 62, pp. 541–552, Jul 1990.
[82] T. Graß, C. Muschik, A. Celi, R. W. Chhajlany, and M. Lewenstein, “Syn-
thetic magnetic fluxes and topological order in one-dimensional spin systems,”
Phys. Rev. A, vol. 91, p. 063612, Jun 2015.
[83] E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl,
P. Hauke, M. Dalmonte, T. Monz, P. Zoller, and R. Blatt, “Real-time dynam-
ics of lattice gauge theories with a few-qubit quantum computer,” Nature,
vol. 534, pp. 516–519, June 2016.
[84] D. V. Else, B. Bauer, and C. Nayak, “Floquet time crystals.” arXiv:1603.08001
[cond-mat.dis-nn], 2016.
[85] T. Lanting, A. J. Przybysz, A. Y. Smirnov, F. M. Spedalieri, M. H. Amin,
A. J. Berkley, R. Harris, F. Altomare, S. Boixo, P. Bunyk, N. Dickson, C. En-
derud, J. P. Hilton, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizin-
sky, R. Neufeld, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva,
S. Uchaikin, A. B. Wilson, and G. Rose, “Entanglement in a quantum anneal-
ing processor,” Phys. Rev. X, vol. 4, p. 021041, May 2014.
[86] T. Albash, W. Vinci, A. Mishra, P. A. Warburton, and D. A. Lidar, “Consis-
tency tests of classical and quantum models for a quantum annealer,” Phys.
Rev. A, vol. 91, p. 042314, Apr 2015.
[87] T. Albash, I. Hen, F. M. Spedalieri, and D. A. Lidar, “Reexamination of the
evidence for entanglement in a quantum annealer,” Phys. Rev. A, vol. 92,
p. 062328, Dec 2015.
[88] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Mar-
tinis, D. A. Lidar, and M. Troyer, “Defining and detecting quantum speedup,”
Science, 2014.
[89] V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy,
J. Martinis, and H. Neven, “What is the computational value of finite range
tunneling?.” arXiv:1512.02206, 2015.
[90] R. Islam, E. E. Edwards, K. Kim, S. Korenblit, C. Noh, H. Carmichael, G.-D.
Lin, L.-M. Duan, C.-C. J. Wang, J. K. Freericks, and C. Monroe, “Onset of
a quantum phase transition with a trapped ion quantum simulator,” Nature
Communications, vol. 2, p. 377, 2011.
[91] P. Jurcevic, P. Hauke, C. Maier, C. Hempel, B. P. Lanyon, R. Blatt, and C. F.
Roos, “Spectroscopy of interacting quasiparticles in trapped ions,” Phys. Rev.
Lett., vol. 115, p. 100501, Sep 2015.
140
[92] H. Haeffner, W. Haensel, C. F. Roos, J. Benhelm, D. C. al kar, M. Chwalla,
T. Koerber, U. D. Rapol, M. Riebe, P. O. Schmidt, C. Becher, O. Guehne,
W. Duer, and R. Blatt, “Scalable multiparticle entanglement of trapped ions,”
Nature, vol. 438, p. 643, 2005.
[93] J.-W. Pan, Z.-B. Chen, C.-Y. Lu, H. Weinfurter, A. Zeilinger, and
M. Z˙ukowski, “Multiphoton entanglement and interferometry,” Rev. Mod.
Phys., vol. 84, pp. 777–838, May 2012.
[94] O. Gu¨hne and G. To´th, “Entanglement detection,” Physics Reports, vol. 474,
no. 1-6, pp. 1 – 75, 2009.
[95] A. M. Kaufman, M. E. Tai., A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss.,
and M. Greiner, “Quantum thermalization through entanglement in an iso-
lated many-body system.” arXiv:1603.04409.
[96] C. Schwemmer, G. To´th, A. Niggebaum, T. Moroder, D. Gross, O. Gu¨hne,
and H. Weinfurter, “Experimental comparison of efficient tomography schemes
for a six-qubit state,” Phys. Rev. Lett., vol. 113, p. 040503, Jul 2014.
[97] C. A. Riofrio, D. Gross, S. T. Flammia, T. M. D. Nigg, R. Blatt, and J. Eis-
ert, “Experimental quantum compressed sensing for a seven-qubit system.”
Preprint.
[98] O. Marty, M. Cramer, and M. B. Plenio, “Practical entanglement estimation
for spin-system quantum simulators,” Phys. Rev. Lett., vol. 116, p. 105301,
Mar 2016.
[99] W.-B. Gao, C.-Y. Lu, X.-C. Yao, P. Xu, O. Guhne, A. Goebel, Y.-A. Chen,
C.-Z. Peng, Z.-B. Chen, and J.-W. Pan, “Experimental demonstration of a
hyper-entangled ten-qubit schrodinger cat state,” Nat Phys, vol. 6, pp. 331–
335, May 2010.
[100] T. Monz, P. Schindler, J. T. Barreiro, M. Chwalla, D. Nigg, W. A. Coish,
M. Harlander, W. Haensel, M. Hennrich, and R. Blatt, “14-qubit entangle-
ment: generation and coherence,” Phys. Rev. Lett., vol. 106, p. 1304, 2011.
[101] D. J. Wineland, J. J. Bollinger, W. M. Itano, and F. L. Moore, “Spin squeezing
and reduced quantum noise in spectroscopy,” Phys. Rev. A, vol. 46, p. R6797,
1992.
[102] M. Kitagawa and M. Ueda, “Squeezed spin states,” Phys. Rev. A, vol. 47,
pp. 5138–5143, Jun 1993.
[103] K. C. Cox, G. P. Greve, J. M. Weiner, and J. K. Thompson, “Deterministic
squeezed states with collective measurements and feedback,” Phys. Rev. Lett.,
vol. 116, p. 093602, Mar 2016.
141
[104] O. Hosten, N. J. Engelsen, R. Krishnakumar, and M. A. Kasevich, “Measure-
ment noise 100 times lower than the quantum-projection limit using entangled
atoms,” Nature, vol. 529, pp. 505–508, Jan. 2016.
[105] A. S. Sørensen and K. Mølmer, “Entanglement and extreme spin squeezing,”
Phys. Rev. Lett., vol. 86, pp. 4431–4434, May 2001.
[106] J. G. Bohnet, B. C. Sawyer, J. W. Britton, M. L. Wall, A. M. Rey, M. Foss-
Feig, and J. J. Bollinger, “Quantum spin dynamics and entanglement gener-
ation with hundreds of trapped ions,” Science, vol. 352, no. 6291, pp. 1297–
1301, 2016.
[107] O. Marty and M. Plenio, “Measuring entanglement in the long-range ising
model.” Private communication, November 2015.
[108] P. Hyllus, W. Laskowski, R. Krischek, C. Schwemmer, W. Wieczorek, H. We-
infurter, L. Pezze´, and A. Smerzi, “Fisher information and multiparticle en-
tanglement,” Phys. Rev. A, vol. 85, p. 022321, Feb 2012.
[109] G. To´th, “Multipartite entanglement and high precision metrology,” Phys.
Rev. A, vol. 85, p. 022322, 2012.
[110] L. Pezze´ and A. Smerzi, “Entanglement, nonlinear dynamics, and the heisen-
berg limit,” Phys. Rev. Lett., vol. 102, p. 100401, Mar 2009.
[111] P. Hyllus, O. Gu¨hne, and A. Smerzi, “Not all pure entangled states are useful
for sub-shot-noise interferometry,” Phys. Rev. A, vol. 82, p. 012337, Jul 2010.
[112] H. Strobel, W. Muessel, D. Linnemann, T. Zibold, D. B. Hume, L. Pezze`,
A. Smerzi, and M. K. Oberthaler, “Fisher information and entanglement of
non-gaussian spin states,” Science, vol. 345, no. 6195, pp. 424–427, 2014.
[113] E. Lieb and D. Robinson, “The finite group velocity of quantum spin systems,”
Commun. Math. Phys., vol. 28, p. 251, 1972.
[114] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalatorre, “Colloquium:
Nonequilibrium dynamics of closed interacting quantum systems,” Rev. Mod.
Phys., vol. 83, p. 863, 2011.
[115] M. A. Cazalilla and M. Rigol, “Focus on dynamics and thermalization in
isolated quantum many-body systems,” New Journal of Physics, vol. 12, no. 5,
p. 055006, 2010.
[116] T. Kinoshita, T. Wenger, and D. S. Weiss, “A quantum newton’s cradle,”
Nature, vol. 440, p. 900, 2006.
[117] S. R. Manmana, S. Wessel, R. M. Noack, and A. Muramatsu, “Strongly corre-
lated fermions after a quantum quench,” Phys. Rev. Lett., vol. 98, p. 210405,
May 2007.
142
[118] M. Babadi, E. Demler, and M. Knap, “Far-from-equilibrium field theory of
many-body quantum spin systems: Prethermalization and relaxation of spin
spiral states in three dimensions,” Phys. Rev. X, vol. 5, p. 041005, Oct 2015.
[119] P. W. Anderson, “Absence of diffusion in certain random lattices,” Phys. Rev.,
vol. 109, p. 1492, 1958.
[120] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, “Possible experimental man-
ifestations of the many-body localization,” Phys. Rev. B, vol. 76, p. 052203,
2007.
[121] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl,
I. Mazets, D. A. Smith, E. Demler, and J. Schmiedmayer, “Relaxation and
prethermalization in an isolated quantum system,” Science, vol. 337, p. 1318,
2012.
[122] D. S. Wiersma, P. Bartolini, A. Lagendijk, and R. Righini, “Localization of
light in a disordered medium,” Nature, vol. 390, p. 671, 1997.
[123] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Cle´ment,
L. Sanchez-Palencia, P. Bouyer, and A. Aspect, “Direct observation of ander-
son localization of matter waves in a controlled disorder,” Nature, vol. 453,
p. 891, 2008.
[124] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Mod-
ugno, M. Modugno, and M. Inguscio, “Anderson localization of a non-
interacting bose-einstein condensate,” Nature, vol. 453, p. 895, 2008.
[125] V. Oganesyan and D. A. Huse, “Localization of interacting fermions at high
temperature,” Phys. Rev. B, vol. 75, p. 155111, 2007.
[126] A. Pal and D. A. Huse, “Many-body localization phase transition,” Phys. Rev.
B, vol. 82, p. 174411, 2010.
[127] M. Serbyn, Z. Papic´, and D. A. Abanin, “Local conservation laws and the
structure of the many-body localized states,” Phys. Rev. Lett., vol. 111,
p. 127201, 2013.
[128] N. Yao, C. Laumann, S. Gopalakrishnan, M. Knap, M. Mu¨ller, E. Demler,
and M. Lukin, “Many-body localization in dipolar systems,” Phys. Rev. Lett.,
vol. 113, p. 243002, 2014.
[129] M. van den Worm, B. Sawyer, J. Bollinger, and M. Kastner, “Relaxation
timescales and decay of correlations in a long-range interacting quantum sim-
ulator,” New J. Phys., vol. 15, p. 083007, 2013.
[130] P. Reimann, “Foundation of statistical mechanics under experimentally real-
istic conditions,” Phys. Rev. Lett., vol. 101, p. 190403, Nov 2008.
143
[131] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, “Relaxation in a com-
pletely integrable many-body quantum system: An ab initio study of the
dynamics of the highly excited states of 1d lattice hard-core bosons,” Phys.
Rev. Lett., vol. 98, p. 050405, 2007.
[132] E. Ilievski, J. D. Nardis, B. Wouters, J.-S. Caux, F. H. Essler, and T. Prosen,
“Complete generalized gibbs ensembles in an interacting theory,” Phys. Rev.
Lett., 2015.
[133] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert,
W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, “Experimen-
tal observation of a generalized gibbs ensemble,” Science, vol. 348, no. 6231,
pp. 207–211, 2015.
[134] J. Berges, S. Borsa´nyi, and C. Wetterich, “Prethermalization,” Phys. Rev.
Lett., vol. 93, p. 142002, Sep 2004.
[135] E. Ising, “Contribution to the theory of ferromagnetism,” Z. Phys., 1925.
[136] M. Srednicki, “Chaos and quantum thermalization,” Phys. Rev. E, vol. 50,
p. 888, 1994.
[137] M. Rigol, V. Dunjko, and M. Olshanii, “Thermalization and its mechanism
for generic isolated quantum systems,” Nature, vol. 452, p. 854, 2008.
[138] G. D. Palma, A. Serafini, V. Giovannetti, and M. Cramer, “Necessity of eigen-
state thermalization,” Physical Review Letters, 2015.
[139] E. H. Lieb and D. W. Robinson, “The finite group velocity of quantum spin
systems,” Communications in Mathematical Physics, 1972.
[140] M. Kollar, F. A. Wolf, and M. Eckstein, “Generalized gibbs ensemble pre-
diction of prethermalization plateaus and their relation to nonthermal steady
states in integrable systems,” Phys. Rev. B, vol. 84, p. 054304, Aug 2011.
[141] E. T. Jaynes, “Information theory and statistical mechanics. ii,” Physical Re-
view, 1957.
[142] M. A. Cazalilla, “Effect of suddenly turning on interactions in the luttinger
model,” Physicals Review Letters, 2006.
[143] J.-S. Caux and R. M. Konik, “Constructing the generalized gibbs ensemble
after a quantum quench,” Physical Review Letters, 2012.
144