Machine-Learning-Based Design of Quantum Systems for
Extreme Sensing
byLiang-Ying Chih
B.S., Physics, National Taiwan University, 2017

A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillmentof the requirements for the degree ofDoctor of Philosophy
Department of Physics
2022

Committee Members:
Murray J. Holland, Chair
Dana Z. Anderson
Graeme Smith
Scott E. Parker
Joshua Combes

iiChih, Liang-Ying (Ph.D., Physics)
Machine-Learning-Based Design of Quantum Systems for Extreme Sensing
Thesis directed by Prof. Murray J. Holland

The constant improvement in metrological devices has advanced the frontiers of science. To
push the limits of precision, the next generation of experiments will be comprised of highly engineered quantum systems. The intrinsic complexity of quantum mechanics often makes engineering
quantum systems to perform a specified task a challenging design problem. We explore a newparadigm for the design and control of devices for performing quantum metrology tasks at levels
that surpass what is achievable by conventional methods. In this thesis, we present three levelsof design philosophy. The first is the periodic modulation of a control parameter that drives thesystem at a resonant frequency, leading to the amplification of a specific signal. We use the generation of momentum-squeezed particle-pairs from a driven Bose-Einstein condensate as an example.
In the second level, we take advantage of a machine learning technique, reinforcement learning, tooptimize the control of quantum states and operations. This approach is demonstrated throughthe design of the matterwave-optic components used for interferometry in a shaken optical lattice.
Going one step further to the third level, we apply reinforcement learning in an end-to-end manner
to design the whole protocol for a metrology purpose. This idea is illustrated through the designof a device that measures rotational signals, i.e., a gyroscope, in a two-dimensional shaken opticallattice. The design is not constrained to the conventional configuration, and leads to a remarkablegain in sensitivity. These design methodologies open new avenues to discovery with the potential
to significantly advance the state-of-the-art quantum sensors.

iii

Acknowledgements

I am lucky to have spent five wonderful years in JILA, where I was provided with a comfortable
research environment and met many awesome people. First of all, I am truly grateful to my advisor,
Murray Holland, for kindly supporting me throughout my PhD and offering unlimited freedom forme to explore different topics. His endless curiosity and creativity has had great influence on me
as a thinker. I am fortunate to have him as a mentor and a friend.
I would like to thank everyone in the Holland group for being fun to work with. I have always
enjoyed our conversations, whether they are about physics or not. In particular, I have learned alot from stimulating discussions with Haonan Liu, Athreya Shankar, John Bartolotta, Simon Jäger,
Jarrod Reilly, and John Wilson. It has also been a pleasure working with Livia Gutierre, Gage
Harmon, and Arnon Goldberg.
I would like to thank everyone on the gyro team, including Josh Combes, Penina Axelrad,
Marco Nicotra, Alex Zouyula, Catie DeLesma, Kendall Mehling, Margaret Raybek, and Jieqiu
Shao for offering a wide range of expertise. Special thanks go to our team leader, Dana Anderson,
for leading interesting and fruitful discussions in the weekly gyro meeting.
I would like to thank friends from Boulder, including Yu-Ju Lee, Xin Xie, Kia Boon Ng,
Yu-Ping Lin, Yang-Zhi Chou, Tianhan Lu, Chih-Ting Hsu, and many more. Not only have they
offered a lot of help, but their company has also made my PhD life more colorful.
I owe a lot of what I have accomplished today to my high school physics teacher, Yi-An Chen.
Without his inspiration and encouragement, I would not have embarked on a journey in physics.
Lastly but not least, I would like to thank my parents for their support in my life.

iv

Contents

Chapter
1 Introduction

1

2 Quantum metrology

5

2.1

Bayes theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.2

Ramsey interferometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.3

Fisher information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.4

Heisenberg limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92.5Squeezed states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.6

Mach-Zehnder interferometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.7

Matterwave interferometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.8

2.7.1

Accelerometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.7.2

Bragg interferometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Gyroscope and Sagnac effect

3 Machine learning. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
17

3.1

Three branches of machine learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2

Reinforcement learning framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3

Tabular Q-learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.4

Deep learning basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4.1

Neural network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

v

3.5

3.4.2

Loss function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.4.3

Gradient descent and backpropagation . . . . . . . . . . . . . . . . . . . . . . 27

Double deep Q-learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.5.1

Improving the learning process . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.6

Variants of reinforcement learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.7

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4 Driving quantum-correlated atom pairs in Bose-Einstein condensates

34

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

4.2

Formalism: General many-body field theory . . . . . . . . . . . . . . . . . . . . . . . 37

4.3

Renormalization of the Scattering Potential . . . . . . . . . . . . . . . . . . . . . . . 40

4.4

Self-consistent ground state solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4.5

Dynamics of the time-dependent HFB system . . . . . . . . . . . . . . . . . . . . . . 48

4.6

Dynamically Generating Squeezed Quasiparticle States . . . . . . . . . . . . . . . . . 52

4.7Interferometry with squeezed quasiparticles . . . . . . . . . . . . . . . . . . . . . . . 56

4.8

Quasi-2D system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.9

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5 Reinforcement-learning-based matter-wave interferometer in a shaken optical lattice

71

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

5.2

Physical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.3

Model-free reinforcement learning for design . . . . . . . . . . . . . . . . . . . . . . . 74

5.4

Component design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.4.1

Beam splitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.4.2

Details of reinforcement learning for target state . . . . . . . . . . . . . . . . 82

5.4.3

Mirror . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

5.4.4

Details of reinforcement learning for target unitary . . . . . . . . . . . . . . . 85

5.4.5

Matter-wave interferometer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

vi
5.5

Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

5.6

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

6 Practical considerations for the shaken lattice accelerometer

95

6.1

From natural units to experimental units. . . . . . . . . . . . . . . . . . . . . . . . 95

6.2

Noise sources in the experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.2.1

Laser linewidth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

6.2.2

Laser intensity fluctuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

6.2.3

Quantum efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

6.2.4

Dark counts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

6.3

Dynamic range . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

6.4

The ultimate limit for sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

6.5

Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

7 How to train your gyro: Reinforcement learning for rotation sensing with a shaken optical
Lattice108

7.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

7.2

Physical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

7.3

Two-path Sagnac matter-wave interferometer . . . . . . . . . . . . . . . . . . . . . . 111

7.4

Effective 1D model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

7.5

End-to-end design with reinforcement learning . . . . . . . . . . . . . . . . . . . . . 118

7.6

Results for end-to-end design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

7.7

Conclusion and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

8 Summary and outlook

125

vii
Bibliography

127

Appendix
A Numerical methods for the Bose-Einstein condensate

134

A.1 Imaginary-time propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
A.2 Spatial derivative in 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
A.3 Spatial derivative in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
B Simulation techniques for the shaken lattice137

B.1 Bloch’s theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
B.2 Momentum basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
B.3 Accelerating reference frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

viii

Tables

Table
5.1

Choice of hyperparameters for the splitting task. The exploration probability ϵdecays linearly for each episode at the rate represented by ϵ decay. The optimizer
we use, ‘Adam’ (which stands for adaptive momentum estimation), is a method for
efficient stochastic gradient descent [42]. . . . . . . . . . . . . . . . . . . . . . . . . . 83

5.2

Choice of hyperparameters for the reflection task. . . . . . . . . . . . . . . . . . . . . 87

ix

FiguresFigure2.1(a) An illustration the Bloch vector (red arrow) for the state cos θ |↑⟩ + eiϕ sin θ |↓⟩ on
the Bloch sphere. (b) An illustration of the Ramsey sequence on the Bloch sphere.
From left to right, the state is initialized in the |↓⟩ state, then a π/2 is appliedthat rotates the Bloch vector onto the equator. The signal couples to the spin asa precession about the z-axis by an angle ϕ, and finally the angle is converted to a
population measurement. (Figure credit: Athreya Shankar.) . . . . . . . . . . . . . .

2.2

7

A sequence that generates the N00N state in a 4-spin system, represented as the
Husimi Q-functions. (a) A coherent spin state pointing in the y-direction. (b)Starting with the coherent spin state shown in (a) and applying the one-axis twisting
Hamiltonian for a duration τ such that χτ = π/2, we get the maximally-squeezed
state. (c) Rotating the spin squeezed state in (b) along the x-axis for angle π/2, we
get the N00N state. The important characteristic of the N00N state is that there are
N fringes lying on the equator, and therefore one can discriminate a rotation angle
that is N times smaller than that for a coherent spin state. This advantage comesat the cost of the reduced dynamic range due to the repetition in the distribution. . 12

2.3(a) The optical Mach-Zehnder interferometer. (b) The Bragg interferometer, whichis the matter-wave analog of the Mach-Zehnder configuration. . . . . . . . . . . . . . 13

x
3.1

The fundamental feedback loop in reinforcement learning. The agent selects anaction based on the current state, and sends the action to the environment, where
the action is executed. The environment then outputs the updated state, along witha reward that indicates how ideal the environmental configuration is after the update. 20

3.2(a) A visualization of the Q-table for the game of tic-tac-toe. (b) The rates of
winning, losing, and tying for the agent that plays first as a function of the total
training episode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4.1

Quantum depletion

1
N

´


GN (x, x) dx as a function of time (the proportion of non-condensate atoms at zero-temperature) simulated with the gapless HFB theory (i.e.,
using Eqs.(4.16), (4.18), (4.21), (4.24), and (4.26)) gives a depletion proportion that
initially scales as ∼ t2 .
4.2

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

Solutions for a system with total atom number N = 6 × 105 in a 1D infinite potential
2 /L.
well of size L, with the scattering potential between atoms given by a = 10−4 l⊥(a) Ground state condensate density found from the gapped self-consistent generalized GPE theory (i.e., replacing Eq. (4.16) with Eq. (4.28) and Eq. (4.18) with
Eq. (4.27)). The length scale over which the condensate density falls to zero at theedges of the box is known is the healing length. (b) Solution to the normal density
GN (x, x′ ) in the ground state as found from the self-consistent HFB theory. (c) Solution to the absolute value of the anomalous density |GA (x, x′ )| in the ground state
as found from the self-consistent HFB theory.

. . . . . . . . . . . . . . . . . . . . . 49

xi
4.3

Snapshots of (a) absolute value of the matrix elements of the upper-left block of
the population matrix, |p|, and (b) absolute value of the matrix elements of theupper-right block of the population matrix, |q|, at t = 0.1 (mL2 /ℏ), starting fromthe initial condition shown in figure 4.2 and then continuously driven with ampli2 /L and frequency ω = 1000 (ℏ/mL2 ). This frequency resonates
tude aac = 10−4 l⊥with the quasiparticles with energy ϵk∗ = 500 (ℏ2 /mL2 ), corresponding to the resonant wavenumbers shown for reference as white lines (at k-index (πk)2 ≈ 500). At
the resonant quasiparticle excitation a clear spike is evident. (c) Density profile ofthe condensate, revealing in general form the spatial dependence of the eigenmodefunction of the resonant quasiparticle excitation. . . . . . . . . . . . . . . . . . . . . 51
4.4

Population in the resonant quasiparticle mode (labelled k) driven with frequency ω =5500 (ℏ/mL2 ) as a function of time for different modulation amplitudes. The amplitudes are Vref (yellow), 2Vref (red), and 4Vref (blue), where Vref = 1.25×10−4 (ℏ2 /mL).
The slope of log pkk at large time is equal to twice the squeezing rate, which is
proportional to the modulation amplitude as shown in Eq. (4.38). The squeezing
rates calculated from the slopes of the curves in the interval t = [0.3, 0.4] (mL2 /ℏ)
are 4.3, 7.2, 14.1 (ℏ/mL2 ), and the squeezing rates calculated from Eq. (4.38) are
3.6, 7.1, 14.3 (ℏ/mL2 ) respectively for amplitudes Vref , 2Vref , 4Vref . . . . . . . . . . 55

4.5

The variance of the quadrature as a function of the angle, at t = 0.05 (mL2 /ℏ) withmodulation amplitude Vac = 2 × 10−3 (ℏ2 /mL), evaluated using Eq. (4.40). The
dotted line is the standard quantum limit, where the variance is equal to 1. For acertain range of angles, the variance falls below the standard quantum limit. Morespecifically, at θ = 0.88π the variance has minimum, which means measurements of
the quadrature along this direction will have the greatest precision. . . . . . . . . . . 57

xii
4.6

The two methods for probing the phase of the quasiparticle squeezing. In both cases,
the scattering length at first oscillates at frequency ω = 5500 (ℏ/mL2 ) for a periodof time τ . At t = τ , in (a) the amplitude of the oscillation remains at zero for atime interval ∆t, and then the scattering length again oscillates for another periodof time τ , and in (c) a phase shift θ is applied to the oscillation . Panels (b) and (d)show the resulting non-condensate fraction at the final time as a function of ∆t orphase shift θ, respectively, for the two cases. . . . . . . . . . . . . . . . . . . . . . . . 60

4.7

Ground state solutions for a system with total atom number N = 6 × 105 in an
infinite circular box of radius R0 . The scattering length is set to a = 3.99 × 10−5 l⊥ .
(a) Condensate density as a function of radial position. (b) Quantum depletion
density as a function of radial position. (c) Anomalous density with the center ofmass position at the center of the trap i.e., r1 + r2 = 0, as a function of the relative
distance |r1 − r2 |. The divergence that scales as the Hankel function of the first kindclose to the origin is a result of 2D scattering theory. . . . . . . . . . . . . . . . . . . 63

4.8

Normal and anomalous densities in a system that is prepared in the ground stateof the self-consistent HFB solution with a positive scattering potential adc = 3.99 ×
10−5 l⊥ , and then subjected to the modulating drive with angular frequency ω =1200 (ℏ/mR02 ) and constant amplitude aac = 3.99 × 10−5 l⊥ . (a) Non-condensate
density as a function of radial potition r = |r1 + r2 |/2 (the origin is on the right)
and time, i.e., GN (r, |r1 − r2 | = 0, t). (b) Magnitude of the anomalous density as
functions of the relative coordinate |r1 −r2 | (the origin is on the right) and time, i.e.,|GA (r = 0, |r1 − r2 |, t)|. Only a small time interval beginning at t = 0.6 (mR02 /ℏ) isshown, so that the oscillations are resolved. . . . . . . . . . . . . . . . . . . . . . . . 64

xiii
4.9

Densities of the condensate and the non-condensate at t = 0.04 and 0.08 (mR02 /ℏ).
The dark blue area is the range of simulation, and the white line indicates thetrapping potential. The yellow arrows represent the density currents, with theirlengths proportional to the magnitude. The system starts with the ground state with
adc = 1.99 × 10−5 l⊥ and then the scattering length is modulated with amplitude
aac = 1.99 × 10−4 l⊥ for all time. The smoothing parameter of the finite circular well
was set to ζ = 0.2. At t = 0.04 (mR02 /ℏ), the atoms in the condensate are pushed
towards the edge of the trap, but because they are of low energy, they cannot escape
the trap. The number of high momentum atoms is still relatively small, and therefore
the density current is also small. The density currents at r = 0.4 R0 point outwards
while those at r = 0.8 R0 point inwards. Outside the trap, there are few atoms,
and the currents are essentially zero. At t = 0.08 (mR02 /ℏ), a large fraction of noncondensate atoms with high energy escape the trap, and the density currents near
the edge of the trap and outside the trap point outwards with large magnitudes. For
clarity, the density currents at t = 0.04 (mR02 /ℏ) are scaled up 3 times compared tothose at t = 0.08 (mR02 /ℏ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

4.10 The size of the non-condensate as a function of time. The size is defined as theradius R that encircles 96.2% of the non-condensate atoms. The system starts withinitial chemical potential µ = 26 (ℏ2 /mR02 ) and is driven with modulation frequency
ω = 1000 (ℏ/mR02 ). The speed of the propagation v is given by the slope of the fittedpline (red), v = 30.8 (ℏ/mR0 ). This value is approximately equal to (ℏω − 2µ)/m,
the speed of a particle whose kinetic energy is half the photon energy minus thechemical potential. We have excluded the data points at times when the motion is
not ballistic, including at early times where few non-condensate atoms escape thetrap and the interaction between the condensate and the non-condensate is significant. 69

xiv
5.1

Interferometers are composed of (i) a beam splitter, (ii) mirrors, and (iii) a recombiner. Examples shown are; (a) an optical Mach-Zehnder interferometer (using
half-silvered or conventional mirrors), (b) a Bragg interferometer (using three shortlight pulses of varying pulse area, π/2 or π), and (c) a shaken lattice interferometer. In (c), the shaken lattice mimics the interferometer components by splitting,reflecting, and recombining the atoms, through design of a specific shaking function
for each case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.2

Framework of reinforcement learning. The agent chooses an action, the environment
responds to the action, and gives the next state and a reward as feedback. In our design task, the action is the translation of the optical lattice, the environment evolves
according to the Schrödinger equation, the observation is the momentum population
distribution, and the reward is a function of the quantum fidelity. Illustrated on theleft is an example neural network as the decision-making agent. The neural network
takes the state represented on its input layer, passes it through a hidden layer, andgenerates a vector of Q-values at its output layer. . . . . . . . . . . . . . . . . . . . . 76

5.3

An example topology of the network that we used in the deep Q-algorithm. . . . . . 79

5.4

Shaking function for splitting. The learned lattice phase, ϕ(t), (top) is allowed totake on one of a discrete set of five possible values that span the range shown. The
momentum probability distribution (bottom) is initialized to the ground state of thelattice and at the terminal time well approximates the desired superposition for abeam splitter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

5.5

Shaking function for reflection. The learned lattice phase, ϕ(t), (top) is sinusoidal
and the amplitude is allowed to take on any one of a discrete set of five values at
each half-cycle. The momentum distribution (bottom) is prepared in −4ℏkL andwell approximates 4ℏkL at the terminal time. . . . . . . . . . . . . . . . . . . . . . . 92

xv
5.6

Time evolution of the matter-wave density throughout the entire interferometry sequence. The white dotted lines separate the plot into five regions. In region I, weapply the splitting protocol. In region II, we allow the matter-wave to propagate
freely in the lattice. The appearance of two wave packets traveling in opposite directions shows that the beam splitter operates as expected. In region III, we apply the
mirror shaking function. The matter-wave undergoes free propagation in region IV
again, and the two wave packets switch directions, demonstrating the functionalityof the mirror. Last, we apply the recombining protocol in region V. Apart fromthe main closed diamond-shaped paths, we observe that there are auxiliary paths
that are fainter but still clearly evident. They arise due to the imperfection of thecomponents, and also due to the side peaks arising from the third-excited Bloch state. 93

5.7(a) Posterior probability density of the acceleration for the first 100 atoms. The true
acceleration that we aim to reveal by the measurements is −3 × 10−3 ωr vr (red line),
and the measurements are sampled from the momentum distribution at the end of
the interferometry sequence, as shown in the inset. (b) Standard deviation of theacceleration estimated using Bayes theorem for up to 104 atoms. We show the resultsfrom both the shaken lattice interferometer (SLI) and the Bragg interferometer andconclude that the shaken lattice interferometer has a higher sensitivity. The standard
deviations σa are roughly inversely proportional to the square-root of the numberof measurements N , for N > 102 . The the black dashed lines are the Cramér-Rao√(CR) lower bounds and scale exactly as 1/ N . . . . . . . . . . . . . . . . . . . . . . 94

6.1(a) Bayesian reconstruction of the true acceleration in the presence of a phase noisethat corresponds to a linewidth ∆ = 2 × 10−4 ωr . The red line indicates the true
acceleration. (b) Standard deviation of the reconstructed acceleration. (c) Meansquared error of the acceleration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

xvi
6.2

(a) Reconstruction of the true acceleration in the presence of intensity noise ξ =2 × 10−4 ωr . (b) Standard deviation of the acceleration. (c) Mean squared error of
the acceleration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

6.3(a) Reconstruction of the acceleration from a system with laser intensity noise as
shown in (c). (b) The standard deviation of the acceleration as calculated from thereconstruction. (c) Time series of the noise. The noise is defined as the ratio between
the intensity fluctuation and the mean intensity. . . . . . . . . . . . . . . . . . . . . 101

6.4(a) Reconstruction of acceleration with E = 0.9. (b) Standard deviation. We assume
that there are originally 10000 atoms, but since E ≠ 1, we only get 9035 counts on
the detector. The inset shows a zoom-in of N = 8000 to N = 9035. . . . . . . . . . . 103

6.5

Error and standard deviation in the presence of dark counts. The absolute value ofthe error grows larger as the noise level increases, whereas the standard deviation
remains almost the same. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

6.6

Momentum distribution for a wide range of acceleration. The free-propagation timeis set to 10 ωr−1 .

6.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104(a) A plot of the estimated acceleration versus the true acceleration. Each acceleration is reconstructed from 10000 samples. The linear relation verifies that theestimations are unbiased. (b) A scatter plot of the error. It can be seen that theerror is uniformly distributed across the whole range of true acceleration. (c) A
histogram of the error, where the width of each bin is 4 × 10−6 (ωr vr ). . . . . . . . . 105

6.8(a) Final momentum distribution for an interferometer where the free propagation
time is 1000 ωr−1 . (b) A zoom-in of (a) near a = 0. (c) The error and the standard
deviation of the estimation of the acceleration from plot (b) with 10000 samples. (d)Histogram of the estimation error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

xvii
7.1

A two-dimensional optical lattice potential in a rotating frame, where atoms aretrapped and their dynamics can be controlled by the shaking of the lattice. The blue
arrows indicate possible trajectories associated with the splitting of the matterwaves. 110

7.2

The clockwise path (blue) and the counter-clockwise path (red) taken by the wave
functions in the two-path shaken lattice gyroscope. The vertical back line is where
the wave functions were initially split and are eventually recombined. . . . . . . . . . 113

7.3

The Fisher information for the machine-learned shaken-lattice gyroscope normalizedwith respect to the Fisher information that results from the gyroscope in which thematterwaves follow the ideal paths of the conventional two-path inteferometer. . . . 115

7.4(a) Average position of the wavepacket in the y-direction. (b) Average momentumof the wavepacket in the y-direction. The lattice depth Vy is set to 10 ωr . . . . . . . 117

7.5

Comparison between the simulations results of the probability distributions in thex direction from the effective 1D model (blue), the full 2D model (orange), and thehighest singular value component from the 2D model (green). . . . . . . . . . . . . . 117

7.6

Phase and the spatial wavefunction as a function of time. The shaking protocol is
generated from the reinforcement learning agent with ϕ̇ as the actions. There are
125 steps in this protocol. The wavefunction is initialized in the ground state witha momentum width σp = 0.1, and evolve with the shaking protocol shown above.

7.7

. 121

Population distribution in the momentum basis and the resulting classical Fisher
information for the machine-learned shaken-lattice gyroscope as a function of therotation rate Ω. The value of the Fisher information is normalized with the Fisher
information of a conventional two-path interferometer with splitting ∆p = 8ℏkL . . . 122

7.8(a) A reconstruction of the rotation rate from 1 sample to 100 samples. (b) The
standard deviation of the rotation rate as a function of the atom number. The resultat large atom number is consistent with the Cramér-Rao bound, which scales as√
1/ N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Chapter 1

Introduction

The constant improvement of ultra-stable lasers with the linewidth now reaching below the
10 mHz level [67], and the ability to cool atoms below the µK regime, have opened up the door
into the quantum engineering realm. The extreme level of control over ultracold quantum systems
has lead to exciting applications such as quantum simulation, quantum computing, and quantum
sensing. In particular, quantum sensing is currently achieving astounding levels of precision andaccuracy, which provides a pathway to challenge the boundaries of physics. For example, we now
have optical clocks that are capable of discriminating the gravitational redshift on the Earth’s
surface at the scale of millimeters [12], and atom interferometers that can support a coherent
superposition of two clouds of atoms even when they are half a meter apart, which allows for thedetection of the gravitational tidal force [47].
One of the key ingredients in cold-atom-based quantum sensors is the Bose-Einstein condensate (BEC). A BEC is the state of matter formed when a cloud of Bosonic particles is cooled
to near absolute zero temperature, such that their de Broglie wavelength is much larger than theinterparticle distance. In this case, their wave functions are mostly overlapped, and can be collectively described by a single macroscopic wave function. Since the matter wave of the BEC ishighly coherent, it is an ideal medium to use for interferometric measurements, in analogy withhow highly coherent lasers are used to carry out precise interferometry.
An example of using BEC to build an interferometer is the experiment by Weidner et al. [100].
In their work, they loaded a BEC into an optical lattice formed by two counter-propagating lasers,

2and sought to construct an accelerometer by ‘shaking’ the lattice with highly engineered patterns.
Through a variety of optimal control methods, they found the shaking functions that gave riseto the splitting, symmetric translation, reflection, and recombination of the matter waves, andcombined these components in the appropriate order for interferometry. The acceleration signals
were then extracted via curve fitting from the interference patterns of the momentum distribution
observed in the time-of-flight images.
When the interactions between atoms in a BEC come into play, even richer and more fascinating dynamics can occur. For instance, in the experiment by Clark et al. [18], they periodicallymodulated the interaction strength between atoms in a BEC, and observed firework patterns formedby the ejected atoms. This experiment is nicknamed the ‘Bose firework’. Intrigued by the experimental observations, we set out to build a model to explain the phenomenon. From our analysis,
we found that the high-momentum atoms generated from the BEC using this relatively straightforward experimental technique are actually entangled in momentum and form metrologically usefulquantum states.
The examples above showcase how high-quality and diverse our toolbox is for the control ofquantum systems. Nevertheless, designing quantum systems for a specific application is still a hard
problem, due to the multiple interference paths or the exponential scaling of states. For all but
the simplest quantum systems, it is in general difficult to figure out what is the best strategy to
combine the tools we have to fully leverage the ‘quantumness’.
But what if we ask the computer to do the engineering for us instead? Machine learning hascelebrated spectacular success over the last decade, due to the exponential increase of computational
power and the recent developments of sophisticated algorithms. Specifically, a branch of machinelearning called reinforcement learning has been shown to be useful when a brute-force search ofstrategies is not possible. It has also demonstrated the ability to learn beyond human intuition
and imagination. The design of quantum sensors fits into this paradigm, and therefore we believe
the marriage of quantum sensing and machine learning can become powerful. Moreover, since
the reinforcement learning approach does not rely on knowledge of the underlying structure of

3
the system that we want to control, we think it is likely that this will open a new paradigm toquantum design, enabling scaling up of the complexity, and allowing the control of the real-world
experimental platforms at JILA and elsewhere.
In this thesis, we present three example designs of quantum systems based on ultracold
atoms for the purpose of high-precision sensing. These examples correspond to three levels ofdesign philosophy, starting from the most intuitive one to one that often generates solutions thatcannot easily be anticipated.
First, we investigate a control protocol that generates quasimomentum pairs from a BoseEinstein condensate by the periodic modulation of the scattering length, as demonstrated in the
Bose firework experiment. This protocol allows the selective amplification of quantum fluctuations
with a predetermined momentum, where the momentum value can be spectroscopically tuned.
Although this control protocol corresponds to the most basic design philosophy that does not
involve any learning, it is still interesting in that the atoms generated form momentum-squeezed
matter waves. These squeezed matter waves have reduced noise in one quadrature and could
potentially be injected into an atom interferometer that aims to achieve a noise level below thestandard quantum limit. The simplicity of this design philosophy also allows us to fully understandthe underlying mechanism and explain the results.
Next, we take advantage of reinforcement learning to design the matterwave-optic components
in a one dimensional lattice. Just like in the original shaken lattice experiment of Weidner, we aim
to find the shaking functions that optimize the quality of the components. With reinforcement
learning, we find the control protocol that leads to a larger momentum splitting of the matterwave. The improvement is partly due to the difference in the exploration space for the shaking
functions. Indeed, our shaking functions bear more dramatic variations than those obtained from
optimal control algorithms. This example demonstrates the potential effectiveness of reinforcement
learning for the control of well-defined quantum states and operations.
Moving forward to the third level, we consider applying reinforcement learning in an endto-end manner to design the entire control protocol in a two-dimensional lattice with the goal of

4achieving the optimal precision for a gyroscopic device. This is the stage where the constraints
set by what physicists have learned from the conventional configurations will be broken, and thefundamental limits to the achievable precision will be challenged.
This thesis is structured as follows. In the next two chapters, we will introduce some background knowledge that will help understanding the material that follows. Chapter 2 will focus on
the basic theory of quantum metrology, which explains why quantum sensors can provide advantage
in sensitivity over classical ones. We will also introduce a couple of text-book standard techniques,
and they will serve as the baseline for comparison with the more complex techniques developed
later. As mentioned previously, we have considered the use of machine learning for the design ofquantum metrological devices. Seeing that machine learning is still a somewhat foreign field formany physicists, it will be necessary to provide a comprehensible guide to it. This will be carriedout in Chapter 3, where we give an overview of machine learning before we dive into reinforcement
learning, lay out the algorithms that we employ, and go through a few pedagogical examples. Afterthe necessary background is introduced, we present our original research. In Chapter 4, we describe
our theoretical investigation of the Bose firework. In Chapter 5, we present the design frameworkand the results for the shaken lattice accelerometer. Then in Chapter 6, we extend our study fromthe previous chapter to real-world applications by taking into consideration the geometry of theexperimental apparatus and the noise sources that may affect the performance of the shaken latticeaccelerometer. In Chapter 7, we present the gyroscope designed by the end-to-end reinforcement
learning approach. Finally, we conclude the thesis in Chapter 8 and provide perspectives for futurework.

Chapter 2

Quantum metrology

Quantum metrology is the study of making measurements with quantum systems to achievea precision that cannot be reached by their classical counterparts. The central limit theorem inclassical statistics tells us that the uncertainty of a parameter scales inversely with the square rootof the number of independent measurements taken. By taking advantage of quantum entanglement,
it is, however, possible to improve upon this classical scaling [29].
In this chapter, we will review the theoretical aspect of fundamental ideas in quantum metrology. Understanding statistics is the first step in extracting useful information from measurement
data, and this is what we will start with. Next, we will introduce a few protocols commonly usedin quantum metrology, including Ramsey interferometry and the ‘squeezing’ of quantum states.
Finally, we will present several experimental techniques that put the abstract protocols mentioned
previously into practice. The ideas introduced in this chapter layout the underlying principles ofquantum metrology on which more complicated extensions or variants can be built.

2.1

Bayes theorem
The Bayesian philosophy is founded on the framework in which there is a fixed probability

model for the data, or interchangeably, measurement outcomes, and from the data one can infer theunderlying parameters of that model. Bayesian parameter estimation usually works when we know
exactly what the forward distribution is, that is, we have a model that predicts the probability
distribution of the measurement outcomes for each possible value that the parameter can take.

6
Assume that one has a prior belief of the parameter ϕ, P (ϕ), where P (. . .) refers to the probability,
and the model describes the probability of getting the data m conditioned on the parameter,
i.e., P (m|ϕ). Then, the probability distribution for ϕ based on the measured data is found fromthe conditional probability as
P (ϕ|m) = N −1 P (m|ϕ)P (ϕ) ,where N =

´

(2.1)

P (m|ϕ′ )P (ϕ′ )dϕ′ is a normalizing factor. If we have multiple measurement outcomes,{m1 , . . . , mN }, then we can apply Eq. (2.1) iteratively, using the updated probability as the prior
to the next measurement, and get
P (ϕ|m1 , . . . , mN ) ∝ P (mN |ϕ) . . . P (m1 |ϕ)P (ϕ) .

(2.2)

The proportionality is resolved by having the posterior probability be normalized to unity.
Typically, the probability distribution becomes narrower as the number of measurements, N ,increases. One finds, in practice, that if each measurement is independent, then the standard√
deviation typically converges as 1/ N due to the central limit theorem.

2.2

Ramsey interferometry
Ramsey interferometry a measurement strategy very often applied in atomic physics [45].

Basically, it is made up of a sequence of operations that goes as follows. Initially, the atom isprepared in |↓⟩ state. A π/2 pulse associated with a rotation about the y-axis is applied androtates the state to the |+x⟩ state. Then, the atom is subject to a rotation around the z-axis by anamount that is associated with the signal that we want to measure. After that, another π/2 pulse
is applied, and we can infer the phase accumulated during the free evolution from a population
measurement in the z basis. The Bloch sphere representation of the sequence is shown in Fig. (2.1).
From a general point of view, the Ramsey sequence is made up of three elements: state
preparation, interaction with the field that bears the signal, formally represented as a rotationinduced by a generator, and finally, a unitary transformation to an appropriate basis followed by

7(a) Bloch sphere(b) Ramsey sequence

Figure 2.1: (a) An illustration the Bloch vector (red arrow) for the state cos θ |↑⟩ + eiϕ sin θ |↓⟩ on
the Bloch sphere. (b) An illustration of the Ramsey sequence on the Bloch sphere. From left toright, the state is initialized in the |↓⟩ state, then a π/2 is applied that rotates the Bloch vector
onto the equator. The signal couples to the spin as a precession about the z-axis by an angle ϕ,and finally the angle is converted to a population measurement. (Figure credit: Athreya Shankar.)

8
a projective measurement. A thriving theme in quantum metrology research is to tailor the initialprobe state and the measurement basis to the signal of interest such that the measurements are
optimal with respect to a certain set of metrics, which is a non-trivial task in a multi-partite or
multi-level system.

2.3

Fisher information
Fisher information is a metric that characterizes how much information one can gain about

the value of a parameter from a measurement sampled from a distribution. The classical Fisher
information is defined as
"

2 #
∂
I (c) (ϕ) = Eϕlog P (m|ϕ)∂ϕ
X∂P (m|ϕ) 2
1=,
P (m|ϕ)
∂ϕ
m

(2.3)

where Eϕ is the expectation value evaluated at the parameter ϕ, m denotes all the possible measurement outcomes, and P (m|ϕ) is the conditional probability of getting a specific outcome m.
The Fisher information may also be calculated fromX  ∂ 2 log P (m|ϕ) P (m|ϕ) ,(2.4)
I (c) (ϕ) =−2∂ϕm
P ∂ 2 P (m|ϕ)
which can be derived from Eq. (2.3) by making use of the identity= 0. With am∂ϕ2
sequence of N measurements, {m1 , m2 , . . . , mN }, the summation should be over all the possible
measurement trajectories, weighted by the corresponding probabilities that can be obtained from
Bayes theorem if independence of the measurements is assumed. The resulting Fisher informationis then(c)
IN (ϕ)

=
=

X

...

X

{mN }

{m1 }

X

X

...

−
−

∂ 2 logQi P (mi |ϕ)
∂ϕ2

i

=

{m1 }

i

X ∂ 2 log P (mi |ϕ)  Y∂ϕ2
X X  ∂ 2 log P (mi |ϕ) =−P (mi |ϕ)
∂ϕ2
{mN }

Y

i

i

P (mi |ϕ)P (mi |ϕ)

{mi }

(c)
N I1 (ϕ) ,

(2.5)

9
where we have used the identityP{mN } . . .

P{mi+1 }

P

{mi−1 }

P

{m1 }

Q

j̸=i P (mj |ϕ) = 1 to get fromthe second line to the third line. What the relation proven above tells us is that the Fisher
information for N independent measurements is N times larger than the Fisher information for asingle measurement.
The inverse of the Fisher information gives the lower bound of the variance of an unbiased
estimator for the parameter, i.e.,
σϕ2 ≥ 1/I (c) (ϕ) .

(2.6)

This lower bound is called the Cramér-Rao bound. With this relation, we can see Fisher information
as a metric for sensitivity, since the sensitivity of a metrological device is characterized by how
small the variance is for the estimated parameter. As can be deduced from Eq. (2.5), the standard√
deviation scales with the number of measurements as 1 N . Using a Bayesian estimator, thestandard deviation approaches the Cramér-Rao bound at a large number of measurements.
The quantum analogue of the Fisher information tells us what is the maximum achievable
information that one can extract from a quantum state. Even with an optimal measurement
basis, the classical Fisher information can only reach the upper bound set by the quantum Fisher
information. With a pure state as the probe state, the quantum Fisher information for the rotation
angle can be calculated from the variance of the generator Ĝ with
I

2.4

(Q)

D E D E 2.
(ϕ) = 4 Ĝ2 − Ĝ

(2.7)

Heisenberg limit
Before we dive into the Heisenberg limit, let us define the standard quantum limit. The

standard quantum limit, despite its name, is what can be achieved with classical measurements.
It says that the standard deviation of a phase distribution for independent measurements scales√
as 1/ N . The Heisenberg limit, on the other hand, is the ultimate limit that can be achievedin a multi-partite system with quantum entanglement. The name is motivated by the Heisenberg
uncertainty principle, from which one would expect a general relationship ∆N ∆ϕ ∼ 1. With10
correlated particles, the standard deviation may therefore approach 1/N , where N is the number of particles in the system. This uncertainty principle scaling can also be proven with thequantum Fisher information. To make it more concrete, the Heisenberg limit can be understoodwith a Ramsey interferometry example. With N correlated atoms, such as the ‘N00N state’,


√
⊗N
⊗N
/ 2, the phase accumulated during free evolution is N times larger,|N00N⟩ ≡ |↑⟩
+ |↓⟩




√
⊗N
⊗N
−iN
ϕ/2
iN
ϕ/2
/ 2. As a result, the phase unceri.e., exp −iŜz ϕ |N00N⟩ = e
|↑⟩
+e
|↓⟩
tainty can be reduced by a factor of N , but at the expense of dynamic range due to the aliasing
effect that arises.
Of course, this is a qualitative description as there is no formal definition of the phase
operator. In fact, due to the periodicity of phase, its standard deviation cannot exceed 2π. For
example, it is clear that a Fock state, which has zero variance in number, violates the so-called
Heisenberg inequality of phase and number. Nevertheless, these ideas are useful in many situations.

2.5

Squeezed states
Squeezed states have reduced variance in one observable at the cost of a larger variance inthe conjugate observable. They become metrologically useful when the direction with reduced
variance aligns with the observable one would like to measure. In quantum optics, squeezed statesof light can be generated in a nonlinear crystal by a physical process referred to as parametric down
conversion. The Hamiltonian that describes the process can be written as
Ĥ = −iℏχ 2â − â† 2
2

(2.8)

where â and â† are the annihilation and creation operators for the photon. Deriving the dynamics
from the Hamiltonian in the Heisenberg picture, we get â(t) = â(0) cosh χt + â† (0) sinh χt, andtherefore the average number is ⟨0| â† (t)â(t) |0⟩ = sinh χt2 . We introduce the two quadratures as
X̂1 = â + â† ,X̂2 = −i(â − â† ) ,

(2.9)

11
then we arrive at the result that the variance in X̂1 exponentially increases while in X̂2 the varianceexponentially decreases. For an initial vacuum coherent state, the variances evolve according to
D E
X̂12 = e2χt ,D

E
X̂22 = e−2χt .

(2.10)In SU(2) (atomic) systems, the analogy to the squeezed state of the optical field systems isthe spin-squeezed state, which has reduced variance in a direction of the collective spin vector. Spinsqueezing can be generated by the one-axis twisting Hamiltonian,Ĥ =where Ŝz = 12(i)i=1 σ̂z

PN

ℏχ 2
Ŝ ,
2 z

(2.11)

is the collective spin in the z-direction. The reduced variance in thecollective spin can be visualized by the Husimi Q-function. The Husimi Q-function represents theoverlap of the state we have with the coherent spin state, |θ, ϕ⟩⟩ ≡ |θ, ϕ⟩⊗N , where |θ, ϕ⟩ is the
single-qubit spin state labelled by its spin direction on the Bloch’s sphere. That is,Q(θ, ϕ) = |⟨⟨θ, ϕ| ψ ⟩ |2

(2.12)

When the one-axis twisting Hamiltonian is applied to a coherent spin state pointing in the ydirection for a time duration π/2χ, we get a maximally squeezed state. If we then rotate the spinsto align with the z-axis, the resulting state is the N00N state, up to a constant phase difference.
The resulting state has the maximum possible Fisher information with respect to a rotation along
the z-axis. The Husimi Q-functions of the coherent spin state, the squeezed state, and the N00N
state are shown in Fig. (2.2).

2.6

Mach-Zehnder interferometer
The Mach-Zehnder configuration is commonly used for measuring the phase difference be-

tween two optical paths through interference. First, the light is coherently split into two paths witha beam splitter, then a mirror reflects each path, and finally the two paths rejoin and a anotherbeam splitter recombines the light to generate the interference. This configuration can be seen as a

12

(a)
Coherent|π/2, π/2⟩⟩

spin

state

(b) Spin squeezed state withχτ = π/2(c) N00N state

Figure 2.2: A sequence that generates the N00N state in a 4-spin system, represented as the Husimi
Q-functions. (a) A coherent spin state pointing in the y-direction. (b) Starting with the coherent
spin state shown in (a) and applying the one-axis twisting Hamiltonian for a duration τ such thatχτ = π/2, we get the maximally-squeezed state. (c) Rotating the spin squeezed state in (b) along
the x-axis for angle π/2, we get the N00N state. The important characteristic of the N00N state isthat there are N fringes lying on the equator, and therefore one can discriminate a rotation angle
that is N times smaller than that for a coherent spin state. This advantage comes at the cost of
the reduced dynamic range due to the repetition in the distribution.

13

Figure 2.3: (a) The optical Mach-Zehnder interferometer. (b) The Bragg interferometer, which is
the matter-wave analog of the Mach-Zehnder configuration.

generalization of the Ramsey sequence, and the use of the mirrors is simply to bring the two paths
back together in real space so that they can interfere with each other. The configuration is shownin Fig. (2.3).

2.7

Matterwave interferometer
The matterwave interferometer takes advantage of the wavelike nature of atoms to carryout interferometry. The de Broglie wavelength of a particle is given by λ = h/p, where p is themomentum of the particle, and h is Planck’s constant. One of the benefits of using matterwaves
over electromagnetic waves is that the wavelength of the matterwave can be much shorter than thatof optical light when the momentum is large, which increases the sensitivity in principle. Particles
are subject to gravitational and inertial forces, and therefore matterwave interferometry is usefulfor gravitational or inertial field sensing.

2.7.1

Accelerometer

The phase accumulated in a Mach-Zehnder matterwave accelerometer can be calculated asfollows. Without loss of generality, we assume that the matterwave is split into a superposition

14√of two momentum eigenstates, |ψ(t = 0)⟩ = (|−p0 ⟩ + |+p0 ⟩)/ 2. We assume that the accelerationon the atom is a constant +a (relative to the reference frame that is accelerating in the negative
direction), and in the first free propagation period t = 0 ∼ T , the state evolves top̂2
exp −i− max̂ T /ℏ |ψ(t = 0)⟩2m= exp −ip̂2 T /2mℏ exp (imax̂T /ℏ) exp −aT 2 [p̂2 , x̂]/4ℏ2 exp ima2 T 3 [x̂, [p̂2 , x̂]]/6ℏ3 |ψ(t = 0)⟩




2 3
= exp ima T /3ℏ exp iaT 2 p0 /2ℏ exp −i(p0 + maT )2 T /2mℏ |p0 + maT ⟩

√


2
2
(2.13)+ exp −iaT p0 /2ℏ exp −i(−p0 + maT ) T /2mℏ |−p0 + maT ⟩ / 2 .
The phase difference between the two momentum states in this period is ∆ϕI = −p0 aT 2 /ℏ. Afterwards, a reflection is applied, switching the signs of ±p0 , i.e.,|p0 + maT ⟩ → |−p0 + maT ⟩ ,|−p0 + maT ⟩ → |p0 + maT ⟩ .

(2.14)In the second free propagation period t = T ∼ 2T , we get ∆ϕII = 3p0 aT 2 /ℏ by following thesame derivation in Eq. (2.13). As a result, the total phase difference accumulated during theinterferometry is ∆ϕtot = 2p0 aT 2 /ℏ. Note that the total phase difference, which is proportional tothe area enclosed in the space-time diagram, is proportional to the momentum splitting, 2p0 andquadratic in the free propagation time.

2.7.2

Bragg interferometer

A crucial part of constructing a Mach-Zehnder matterwave interferometer is the capabilities
to coherently split the wavefunction into an equal superposition of two momentum states and swap
the two momenta while preserving the relative phase between the states. The simplest way to do
this is through Bragg diffraction [56], where a standing-wave laser pulse is applied to the atom,and the atom absorbs a photon in one direction, gets excited to a virtual state, and then emits aphoton in the opposite direction and returns to the ground state. This process gives the atom amomentum kick of 2ℏkL , where kL is the wavenumber of the laser.

15
The picture described above works in the Bragg regime, where the laser intensity is notstrong enough to induce higher-order diffraction, i.e., diffraction processes that are not exactly
resonant. With the appropriate choice of reference frame, we can assume that the only relevant
momenta are +ℏkL and −ℏkL , since they have the same kinetic energy and therefore their energies
are degenerate. The Rabi frequency of the laser pulse, Ω, should be smaller than the frequency
corresponding to the energy cap between the |ℏkL ⟩ and the |3ℏkL ⟩ states, i.e., Ω ≪ 4ℏkL2 /m, sothat this transition and the ones above are off the energy shell. These assumptions simplify thefull Hamiltonian to allow an effective two-level Hamiltonian,
ℏΩ
Ĥ =
2




|−ℏkL ⟩ ⟨+ℏkL | + |+ℏkL ⟩ ⟨−ℏkL |.

(2.15)

If we see the two-momentum model as a spin-1/2 system, then this Hamiltonian is essentially
proportional to the Pauli matrix σ̂x . Applying the laser pulse for a duration τ is equivalent torotating the pseudo-spin around the x-axis for an angle Ωτ on the Bloch sphere. This is reminiscentof the first field in Ramsey interferometry as introduced earlier. It is then straightforward to seethat to split the atoms from the initial |−ℏkL ⟩ momentum state to an equal superposition of
the two momentum states, one should apply the laser pulse with integrated intensity such that´τ′′0 Ω(t )dt = π/2. This transfers the pseudo-spin from the north pole to the equator of the
Bloch sphere. Similarly, recombination is achieved by applying another π/2-pulse. To reflect the
momentum states, which is equivalent to swapping the |±ℏkL ⟩ states, one should apply the laser´τpulse such that 0 Ω(t′ )dt′ = π. An illustration of the pulse sequence is shown in Fig. (2.3).
At the beginning of the Bragg interferometry sequence, the momentum of the atom cloud
may not correspond precisely to the desired −ℏkL . If this is the case, one should adjust the laser
frequencies so that the standing wave that results from the interference is moving at the velocity
that matches the center-of-mass velocity after splitting.

16

2.8

Gyroscope and Sagnac effectGyroscopes usually operate based on the Sagnac effect, which gives rise to a phase difference

between two counter-propagating paths when the system is rotating. For matter waves in freespace, the Hamiltonian is
Ĥ =

p̂2
− Ω · L̂ ,2m

(2.16)

where L̂ = r̂ × p̂. Since the first term in the Hamiltonian commutes with the second term, we only
need to consider the phase shift induced by the rotation term. If we adopt the semiclassical picture
and make the substitution p = mṙ, then the phase difference between the two paths that follow
the same trajectory but in the opposite directions is
∆ϕ =
=
=

ˆ
ˆ
1
1
(+)
(+)
mΩ · r× ṙ dt −mΩ · r (−) × ṙ (−) dtℏℏ˛2mΩ · rdrℏ
4mΩ · A.
ℏ

(2.17)In the equation above, (+) and (−) labels the counterclockwise and clockwise paths respectively,
and the contour integral is along the counterclockwise path, which encloses an area A.
Similarly, the Sagnac effect also gives rise to a phase difference between two counter-propagating
paths of optical light in a fiber-optic gyroscope. The phase difference can be deduced from Eq. (2.17)
by seeing the factor m/ℏ as p/ℏv, which is kL /c for light. The result is ∆ϕ = 4kLcΩ·A , where kL is
the laser wavenumber and c is the speed of light.

Chapter 3

Machine learning

Machine learning is the field of study that aims to ‘teach’ machines to make predictions
and decisions based on data and observations without explicitly programming them to do so. Asa powerful tool, machine learning has recently gained more and more attention in the physicalscience communities. One specific branch of machine learning, reinforcement learning, has grownin importance over the last few years with the introduction of a computer program, AlphaZero,which demonstrated the capability to beat the best players in highly non-trivial games like go andchess [85, 87, 86]. The main challenge of developing computer programs that play these games
is that the number of possible board configurations grows exponentially, and it quickly becomes
impossible to determine what moves are the best through an exhausted search. The design andcontrol of physical systems is very much like a game in some sense. The set of metrics associatedwith achieving certain design goals is like the score one gets in a game, and the design with thehighest score is the winner. The games in this context are also highly non-trivial in that there
are usually multiple knobs that can be turned, and the number of combinations for the setup andparameters to be tuned in a time sequence is intractably large. In this chapter, we will briefly
introduce different branches of machine learning and their applications, and then dive deeper intoa specific branch, reinforcement learning. We will start with an the general framework and thendescribe in detail how a specific type of reinforcement learning, Q-learning, works. We will alsocomment on how reinforcement learning can benefit the design of physical systems. The purposeof this chapter is not to present new research in machine learning, but to cover the fundamental

18knowledge one needs to learn to apply machine learning to the development of new designs forquantum systems.

3.1

Three branches of machine learning
Supervised learning is a branch of machine learning that learns from labelled data. One classof problems that supervised learning is used for is regression. An example is learning from historical
data of temperature from a certain region, and trying to predict the temperature in that region byweighing in factors like pressure maps, humidity profiles, and terrain. Another class of problems
is classification, where the machine learns from pre-labelled data to sort new data into categories.
An example is training the machine to recognize hand-written digits with labelled images.
Unsupervised learning refers to a wide variety of machine learning algorithms that learn
from unlabelled data. The goal is to find patterns from the distribution of data. One class ofsuch a problem is clustering. Imagine information from customers that includes age, income, etc.
represented by data points in a high dimensional space. Algorithms like k-means clustering put
similar customers into the same group to get insight as to what products they may like according totheir common traits with other people. In this case, there are no existing samples labelled by whatgroup they belong to. Instead, how they are grouped together depends on the overall distributionof the customer base.
Another class of problem that unsupervised learning deals with is dimension reduction. To
make an example easy to visualize, data points in a three-dimensional space may only be spreadout on a two-dimensional plane and have little variance in the direction perpendicular to the plane.
One can therefore find the two vectors that span the plane by finding the axes with the highest
variances. This method is called the principle component analysis . More complex methods include
the autoencoder, which is a systematic approach to find the latent space representation in the input
data.
Reinforcement learning is a unique branch of machine learning that learns from trials anderrors instead of from pre-existing examples, whether labelled or un-labelled. This is the branch

19that we will focus on in this thesis, and will be elaborated on more in the next section.

3.2

Reinforcement learning framework
Typically, reinforcement learning involves a feedback loop between an agent that invokes

actions based on the observed state of an environment, and an environment that provides rewards
to the agent based on the resulting outcomes. This framework can be formalized as a Markov
decision process, where the probability of transition to the next state is only dependent on thecurrent state and the current action, and not on prior history. The agent is tasked with the goal
to discover the sequence of actions for which it receives the highest possible long-term reward [90](see Fig. 3.1). It does this by trial and error, iteratively improving its actions in such a way as tocorral the environment toward a target configuration.
Evolution in this framework is referred to as a trajectory, where the state is initially prepared,
and then a sequence of actions and corresponding state updates are performed. The trajectory steps
continue until the chain comes to an end when a predetermined condition is reached, which might
be a certain number of steps, or a terminating state of the environment. An important concept
is the idea of an optimum trajectory starting from any state, which is a trajectory such that thelong-term reward is the maximum possible.
Decisions are made according to a policy framework, which we label as π, that instructs theagent how to carry out specific actions based solely on the current environmental state. That is, πtakes an environmental state s and maps it to an action a either deterministically or probabilistically.
The learning paradigm is associated with evolving towards an optimal π that maximizes areturn according to the following algorithm. At each time step, t, the agent receives a reward,
Rt , that is a pure function of the environment and measures the degree to which the environmentis in the desired configuration. We have been using the term ‘long-term reward’, but we refer tois formally defined as the return. The return is the cumulative discounted reward, which can be20

Figure 3.1: The fundamental feedback loop in reinforcement learning. The agent selects an action
based on the current state, and sends the action to the environment, where the action is executed.
The environment then outputs the updated state, along with a reward that indicates how ideal theenvironmental configuration is after the update.

21
expressed in terms of the immediate reward and the future rewards,
∞
X

Gt =

γ k Rt+k .

(3.1)

k=0

This expression accounts for the fact that actions that are taken at time t affect all later times,and γ ∈ [0, 1] is the discount factor that characterizes how much the future rewards are valued dueto the diminishing influence the current action has as the time goes further into the future.
To evaluate the policies, we define a value function, which calculates the expected return
starting from a state s, i.e.,V π (s) = Eπ"∞X#
γ k Rt+k St = s ,

(3.2)

k=0

where the expectation value is over all possible trajectories denoted by the sequential map of statesπ

to actions {St+k −→ At+k }, for k ≥ 0. The optimal value function is the maximum of the valuefunction over all possible policies, i.e.
V ∗ (s) = max V π (s) .
π

(3.3)

Since the policy is an instruction set as to how to select an action, we also need a way toevaluate the desirability of each action. This leads us to defining the action-value function, also
known as the Q-function, where Q stands for quality. It calculates the expected return startingfrom a state s, taking the action a, and then following the policy π, i.e.,"∞#X
Qπ (s, a) = Eπ
γ k Rt+k St = s, At = a .

(3.4)

k=0

A simple relation between the value function and the action-value function isV π (s) =X

π(a|s)Qπ (s, a) .

(3.5)

a∈A

The optimal action-value function isQ∗ (s, a) = max Qπ (s) ,
π

(3.6)

22
and the optimal policy π ∗ can be found by maximizing Q(s, a) over all actions, i.e.
π ∗ (s) = argmax Q(s, a) ,a∗Q= max Qπ (s, a) .

(3.7)

π

By isolating the immediate reward from the future rewards in Eq.(3.2), we can derive the Bellman
expectation equation for the value function,"#∞X
γ k Rt+k St = s
V π (s) = Eπ Rt +


k=1
π


= Eπ Rt + γV (St+1 ) St = s
"
#
X
X
′
π ′
=
π(a|s) r(s, a) + γ
P (s |s, a)V (s ) ,

(3.8)

s′ ∈S

a∈A

where s′ represents the state that follows after the action a is applied to the state s. Similarly, theBellman expectation equation for the action-value function can be written as
Qπ (s, a) = r(s, a) + γEπ [Qπ (s′ , a′ )] .

(3.9)

By recognizing that when the policy is optimum the expected return should be maximized, we
arrive atQ∗ (s, a) = r(s, a) + γ max Eπ [Qπ (s′ , a′ )]π∗= r(s, a) + Es′ [Qπ (s′ , π ∗ (s′ ))]∗= r(s, a) + Es′ [maxQπ (s′ , a′ )] .
′a(3.10)In the scenario where the transition is deterministic, the above equation boils down to the Bellman
optimality equation [7]Q∗ (s, a) = Rt + γ max
Q∗ (s′ , a′ ) ,′a(3.11)

which we use as our fundamental design equation. The learning algorithm we will focus on is Qlearning, which relies on learning to approximate the actual optimal Q-function as a way to obtain
the optimal policy.

23

3.3

Tabular Q-learning
As mentioned previously, learning a good representation of the Q-function is the key to Q-

learning. One of the simplest ways one can think of is to store the Q-value for each state-action pair,Q(s, a), in a lookup table where the rows correspond to the states, s, and the columns correspondto the actions, a. The Q-values in the table are updated through the temporal-difference method,


′
′ ′
Q(s, a) ← Q(s, a) + α r(s ) + γ max
Q(s , a ) − Q(s, a)
a′


′
′ ′
= (1 − α) Q(s, a) + α r(s ) + γ max
Q(s , a ) ,′a(3.12)

where Q(s, a) moves toward the target return, Y (s) ≡ r(s′ ) + γ maxa′ Q(s′ , a′ ) with learning rate α.
As introduced above, machine learning algorithms typically involve many hyperparameters that canbe tuned, making the learning tasks non-trivial. The simplicity of tabular Q-learning allows us tolook under the hood and thereby to gain access to how the exploration for the optimal solution
works. It also has the advantage of being quick, which can become crucial if one were to implement
the learning algorithm as a real-time control method in an actual experiment. Tabular Q-learning
is suitable for small-sized problems, where an exhaustive list of all the possible states can be stored
and processed on the computer with reasonable memory requirements.
Here we show how tabular Q-learning works with a simple example, where we train the computer program to learn to play the game tic-tac-toe. In the program, there are two reinforcementlearning agents that represent players on the two sides, and an environment representing the 3 × 3
board configuration. The values in the Q-tables for both agents are all initialized to 0.5, and theboard in the environment is initialized to be empty. The agents choose the action associated withthe highest value among the values in the row corresponding to the current state in the Q-table.
Every time agent 1 makes a move, placing ‘O’ on the board, the environment evaluates whether theboard configuration is terminal. If not, agent 2 will also make a move, placing ‘X’ on the board,and the resulting board configuration will be sent to this agent 1 as the next state. If the board
configuration is terminal, then a reward is given to the agent. The reward is 1 if the agent wins, 0

24
if the agent loses, and 0.5 if the two agents draw. Once a game is finished, each agent updates theirown Q-table based on the sequence of states and actions stored in their memory. The updates startfrom the last state before the terminal state and the action chosen in the last step. The Q-value
for this state-action value is updated to the reward that the agent received for the terminal state,i.e., Q(sT −1 , aT −1 ) ← r(sT ). Then, going backward in the sequence, the Q-values associated withthe state-action pair at each step is updated following
Q(st , at ) ← (1 − α)Q(st , at ) + αγ max
Q(st+1 , a′ ) .
′a(3.13)

This process is repeated for many iterations until the Q-values in the table converge. We show
the winning rate of agent 1, the winning rate of agent 2, and the rate that the two agents draw in
Fig. (3.2. Initially agent 1 has a higher winning rate because the game favors the player who plays
first if both sides are playing randomly, and after the agents are both fully trained they draw inalmost every game.
Of course, teaching the computer to play tic-tac-toe may not seem all that impressive, as
humans can learn to play this game optimally with few trials. However, a point worth emphasizingis that the agent does not need to know what the environment actually is. To the agent, the actions
are just integers that enumerate elements from a set of possible actions. They are only meaningful
to the agent through the updated state that the environment feeds back. The complete conceptual
separation between the agent and the environment implies that one can potentially solve completely
different problems, which are encoded in the environment, with the same learning agent. This is
the essence of model-free learning. For example, one can replace the environment with a simulationof a physical system, and have the same agent that learned tic-tac-toe learn to control the system toachieve the desired configuration. In fact, an important benefit of having the model-free structure
for physicists is that the environment can also potentially be replaced by an actual experimental
apparatus, which may not be perfectly described by a theoretical simulation or model, but the agent
would still be able to learn from observations from the experiment nonetheless. This is an extremely
powerful concept that is not typically shared by gradient-based optimal control approaches that

25

(a)

(b)

Figure 3.2: (a) A visualization of the Q-table for the game of tic-tac-toe. (b) The rates of winning,losing, and tying for the agent that plays first as a function of the total training episode.

26
require the model to be known.
Understanding the simple tabular Q-learning algorithm is a prerequisite to moving forward
to more complicated agents, such as the double deep Q-learning algorithm that follows.

3.4

Deep learning basics
Deep learning utilizes the cascade of multiple consecutive layers of non-linear functions toapproximate functions with complicated behaviors. In reinforcement learning, the Q-function is
the function we would like to approximate, and deep learning could therefore be useful for thispurpose. In this section, we introduce the key ideas of deep learning, including the neural network
approach, the loss function, and the learning stage through gradient descent.

3.4.1

Neural network

A typical neural network structure is made up of an input layer, several hidden layers, andan output layer. Each layer is connected to the next one by a linear map followed by a non-linear
activation function. This network topology was motivated by the highly connected active neurons
in the brain. The most naı̈ve neural network is the simple feed-forward neural network, where eachnode in a layer is connected to all the nodes in the next layer. In general, the formula for forward
propagation can be written asz (l) = W (l) a(l−1) + b(l)a(l) = g (l) (z (l) ) ,(3.14)

where l labels the layer. In the formula above, W is the weight matrix with dimensionality n × m, b
is a m×1 bias vector, and g(x) represents the non-linear activation function. The weights and biases
are the parameters to be trained in the neural network, and are often denoted by θ collectively.
The initial condition a(0) is given by the sample input vector xi , and a(L) , where L is the total
number of layers, is equivalent to the output vector. Two common activation functions include
ReLU, ReLU(x) = max(x, 0), and the sigmoid function, σ(x) = 1 +1e−x , although there are many

27others used. The use of the non-linear activation function is essential, for if it is omitted, the wholefunction would be a linear transformation, and the multi-layer structure would be redundant. The
neural network forms a highly parameterized hypothesis function, denoted by hθ (x), which we use
to attempt to approximate any actual function with high fidelity.

3.4.2

Loss function

To evaluate the efficacy of the hypothesis, we define a loss function that describes how faithful
the hypothesis is to the desired actual function. The smaller the loss, the better the hypothesis is.
Of course, we do not know the actual function, for that is exactly why learning is carried out in thefirst place. Therefore, we use samples generated from the actual functions and evaluate the loss by
averaging over these samples. One of the most common loss functions is the mean squared error,
L(θ) =1 X|hθ (xi ) − yi |2 ,N

(3.15)iwhere i labels the training sample index, N is the number of samples used to calculate the loss, xi
and yi are the input vector and the target vector respectively for each sample.
Typically, the learning cycle in machine learning algorithms involves techniques to minimize
a loss function by variation of the adjustable parameters θ in order to arrive at the optimal map
for the purpose at hand.

3.4.3

Gradient descent and backpropagation

To minimize the loss function, the most common numerical approach is gradient descent.
With each step of gradient descent, the parameters in the neural network move in the direction of
the locally steepest descent of the loss function. The update rule corresponding to this approach isθ ←θ−α∂L,∂θ(3.16)

where α is the learning rate that determines the step size. If the learning rate is too small, it willtake a long time to converge to the minimum loss. However, if the step size is too large, then thegradient descent process may not even converge at all.

28
Since the neural network has multiple connected layers, the gradients need to be calculatedusing the chain rule of differentiation, and the calculation of the gradients starting from the output
layer all the way to the input layer is called backpropagation.

3.5

Double deep Q-learning
The objective of deep Q-learning is to train the deep Q-network so that it represents thereturn defined in Eq. (3.1) well and therefore satisfies the Bellman optimality equation,
Q(s, a) = r(s′ ) + γ max
Q(s′ , a′ ) .
′a(3.17)

We train the deep Q-network using the double deep Q-learning algorithm [94]. In this algorithm, an extra network called the target network is included to estimate the future return.
The reason to do it this way is that Eq. (3.17) involves a self-consistency element. If during thetraining stage, the Q-values on both sides of the equation are determined by the same Q-network,
then updating the Q-network will simultaneously change both the Q-value Q(s, a) and the return
Y (s′ ) that the Q-value is supposed to converge to. This creates a feedback cycle that can makethe learning process extremely unstable, and in fact the Q-values may never converge. In orderto mitigate this problem, we employ a Q-network, Qθ , just as described, from which the optimal
action is chosen by the agent through their policy. However, we then calculate the expected returnfor the optimization of the Q-network using a Q-value corresponding to this action but obtained
from a second neural network, the target network, Qθ′ . The target network is only updated more
slowly through a weighted average. This avoids the unstable feedback since the Q-network andtarget network are weakly correlated at early stages of the learning process [61]. Furthermore, thismethod resolves issues with overoptimism in the standard deep Q-learning algorithm [94]. A sketchof the algorithm is presented in Algorithm 1.
The algorithm is described as follows. First of all, we initialize the weights and biases in the
Q-network randomly, and initialize the target network such that its network parameters are thesame as those of the Q-network. Then, we iterate over a large number of episodes. In each episode,29

Algorithm 1: Double Deep Q-Learning
Initialize Q-network Qθ and target network Qθ′ ← Qθ ;for trajectory=1:episodes dos ← s0 ;while d=0 do
if rand() > ϵ thena ← arg maxa Qθ (s, a);elsea ← rand(A);end[s′ , r, d] ← environment(s, a);Store (s, a, s′ , r, d) in replay buffer;
Sample from replay buffer: {(si , ai , s′i , ri , di )| i = 1, ..., B};P
L = B1 i [Qθ (si , ai ) − Yi ]2 , Yi = ri + γ Qθ′ (s′i , arg maxa′ Qθ (s′i , a′i ));Update Qθ : θ ← θ − α∇θ L;
Update Qθ′ : θ ′ ← τ θ ′ + (1 − τ )θ;s ← s′ ;
end
end

30
a trajectory of consecutive states and actions is generated. A trajectory starts with an initial states0 and ends when we reach a terminal state or a threshold number of steps.
In each step, we choose the action that maximizes the Q-value for the current state most of
the time, but with probability ϵ we choose a random action. This policy is known as the ϵ-greedy
policy. In practice, ϵ decays from 1 to a number that asymptotically approaches 0 over time, since
exploitation becomes more important than exploration as the learning process proceeds. The actionis then fed into the environment, and the environment outputs the subsequent state s′ , a reward r,
and a Boolean d that indicates whether the trajectory terminates at this step. A tuple (s, a, s′ , r, d),which represents the experience we get from taking this step, is pushed into a memory stack known
as the replay buffer.
To minimize the error of Qθ (s, a) for representing the return, we calculate a loss function
Pdefined as the squared difference between Qθ (s, a) and Y (s′ ) ≡ k=0 γ k R(sk+1 ), where sk represents the state we get after k steps from s . The future return in Y (s′ ) (with k > 0) is evaluatedby the output value of the target network for the next optimal action. The next optimal action
is determined by the action that maximizes the output values of the Q-network, given the next
state s′ as input. Since we need the Q-values to do well for all possible pairs of (s, a), we draw B(the batch size) samples from the replay buffer, where each sample is labelled by an index i, andcalculate their average loss.
The parameters in the Q-network are updated through gradient descent in the direction of
decreasing average loss. The parameters in the target network are updated by taking the weighted
average of the previous parameters from the target network and the updated parameters from the
Q-network.

3.5.1

Improving the learning process

Although we have outlined the learning process in the previous section, there are a few usefultricks to improve the performance of the learning process that are not discussed much in literature.
The first trick is to iterate the process multiple times, each time adding more complexity

31
to it. For example, one can consider dividing the time into more steps. Another example is to
add more actions to the action set. In tabular Q-learning, this simply means that we can growthe Q-table horizontally, adding more columns associated with the new actions to the table. In
deep Q-learning, this involves adding more nodes on the output layer, and initializing the weights
and biases associated with these extra nodes randomly while copying the parameters in the smaller
model trained in the the previous iteration to the new bigger model. The search space for the
optimal protocol is exponentially huge, and one might find learning with a large action set and longsequence difficult. The iterative approach serves as a systematic way to guide the learning process,and often leads to better convergence.
The second trick is inspired by transfer learning, where a model trained for one scenario is
applied to another scenario with some retraining. As an example, one can imagine that a deep neuralnetwork agent is trained first using the computer simulated environment, and then made to interact
with an actual experiment. It is anticipated that the learned protocol would not be optimal for theexperiment since there are usually imperfections in a real-world scenario. However, the simulated
environment should still be close to the experimental environment, and therefore, we could allow
only the last few layers of the neural network to vary and be retrained while fixing the rest of theparameters as static. Since the neural network is close to being fully trained, the retraining process
is expected to take many fewer episodes. Due to the long run time of the experiment comparedto the computer simulation, implementing machine learning on an actual experiment has naturally
been a more difficult task than for the numerical counterpart. This trick could potentially savethe training time on the experiment, and hopefully make reinforcement learning effective in thereal-world scenario.

3.6

Variants of reinforcement learning
We have focused on Q-learning for the most part of this chapter. There are in fact a wide

variety of alternative learning algorithms for the agent, which we will introduce briefly in thissection.

32
In the policy gradient algorithm, the agent makes decisions based on the policy network,
which takes the state s as the input and outputs the probability to select each action a, i.e., π(a|s).
The agent selects a random action at each step with weights based on this probability distribution.
The policy is parametrized by a set of variables denoted by θ, and the learning process entails
maximizing an objective function with respect to θ. A common choice of the objective function isdefined as the expected value function for the initial state, i.e.,J(θ) ≡ Eπθ [V (s0 )] .

(3.18)The expectation value is evaluated over {τ }, a set of trajectories sampled over πθ . Thus, theobjective function can be written asˆJ(θ) =≈πθ (τ )R(τ )dτ
X Xi∈{τ }(i)Rt .

(3.19)t

To maximize the objective function, we calculate the gradient of it with respect to θ,ˆ∇θ J(θ) =≈πθ (τ )∇θ log πθ (τ )R(τ )dτ
X Xi∈{τ }

t

(i)

(i)

∇θ log πθ (at |st )

X

(i)

Rt ,

(3.20)tand update the parameters θ by
θ ← θ + α∇θ J(θ) .

(3.21)

The actor-critic model for the agent can be seen as a hybrid of Q-learning and policy gradient.
In this framework, there is an actor that selects an action based on its policy network, which outputs
the probability distribution of the actions, and a critic that tells the actor how good the selectedaction is based on evaluation of the corresponding Q-value from its Q-network.

3.7

Summary
Machine learning has been demonstrated to be a powerful tool. The power comes at the costof complexity; there are many moving parts in a machine learning algorithm that can be improved,

33
which is also why machine learning is such an active research field. From a broader perspective,
reinforcement learning is the branch of machine learning most relevant to quantum design problems
for complex systems since it learns a pathway through a multi-parameter landscape by trial anderror. In such a situation, a brute-force search would be prohibitive.

Chapter 4

Driving quantum-correlated atom pairs in Bose-Einstein condensates

4.1

Introduction
The ability to tune the two-body scattering length in a Bose-Einstein Condensate (BEC) by

varying the magnitude of a magnetic field in the vicinity of a Feshbach resonance has been employed
in a number of seminal experiments that aim to investigate controlled non-equilibrium quantum
dynamics. One such example is the so-called ‘Bosenova’ experiment by Donley et al. [21], in whicha BEC was subject to a sudden change of the scattering length from a small positive value to a large
negative value. This resulted in a change of the sign of the mean-field interactions from repulsive,
where the gas is mechanically stable, to attractive, where the compressibility may be negative andthe gas is then unstable [79]. What was observed experimentally after this abrupt change in the
scattering length was a collapse and subsequent explosion of the quantum gas in a manner that
resembled an astrophysical supernova. Theoretical models were developed and illustrated that the
emergence of a pairing field in the underlying quantum many-body system can explain the observed
burst of non-condensate atoms [59].
More recently, Bose ‘firework’ experiments [18] have observed pairs of high momentum atoms
emitted as jets from a condensate driven by a periodic modulation of the two-body s-wave scattering
length. These experiments demonstrated a protocol for resonantly amplifying quantum fluctuations
with well-controlled momenta when starting from a stationary BEC. The fact that the jets wereobserved to be correlated in emission direction motivates us to consider whether the many-bodypairing field played an important role in the dynamics, in a similar manner to the Bosenova system

35
previously studied. This aspect is related to other calculations that explore the second ordercoherence of the gas [4].
In the case of a dilute quantum gas, the Gross-Pitaevskii equation (GPE) provides an accurate
description of the equilibrium and time-dependent behavior of the BEC. In this framework, theinteracting condensate is completely described by a mean-field superfluid order parameter. The
GPE framework has been extensively applied to model the behavior of BECs at zero temperature,
and also to their coherent manipulation through externally applied potentials. However, when
there is a significant portion of non-condensate atoms, the GPE will fail to provide an accurate
description of the system. A small amount of non-condensate atoms is always present even atzero temperature in a dilute quantum gas arising from the beyond mean-field fluctuations thatare due to the finite interaction strength. It is of course possible to generate substantial fractionsof non-condensate atoms by driving a pure condensate in a variety of ways, and this is typicallyunavoidable when the currents that generate the magnetic confinement fields contain stochastic
noise. Furthermore, non-condensate atoms are always present in systems with finite temperaturesince they embody the thermal excitations. In order to capture the essential dynamics associatedwith the non-condensed component, a theory that goes beyond the mean-field approximation isnecessary.
A systematic extension to the simplest mean-field approach given by the GPE is the HartreeFock Bogoliubov (HFB) formalism that takes into account the interactions between three components. These components are the condensate, the non-condensate, and the pairing field of the fluctuations. In this formalism, the elementary excitations are described as Bogoliubov quasiparticles[10], and the ground state condensate is a vacuum of such quasiparticles. The quasiparticle creation
operator is a linear combination of the atom creation operator (particles) and the atom annihilation operator (holes). The vacuum state due to interactions possesses a portion of non-condensate
atoms referred to as quantum depletion. The ‘firework’ experiments that tune the scattering length
by applying an appropriate external magnetic field potentially allow all of these components—themean-field, non-condensate, and pairing field—to be controlled, manipulated, and engineered. In

36this chapter, we derive the solutions of the HFB theory as applied to well-controlled experimental
geometries in order to determine the efficacy of this framework for providing a theoretical basis forthe observations from the ‘firework’ experiments.
One approach for describing collective excitations of a condensate involves solving the Bogoliubov de Gennes equations [109], which is most accurate when the excitations are weak. When the
excitations are not weak, and the non-condensate fraction can be significant, a more complete approach must be used such as the time-dependent self-consistent HFB equations [16], and that is themethod we will focus on in this chapter. Note that one alternative approach that can incorporate
the excitations is through the addition of a noise source to introduce fluctuations directly into theGross-Pitaevskii equation [25]. However, this assumes by construction that the many-body state
can be accurately described by a unique macroscopic wavefunction, and therefore a more complete
theory is needed to describe two-body correlations.
We emphasize the importance of the pairing field in our analysis. Pairing gives rise to ananomalous density that allows us to investigate the coherence of the system and explore methods
to probe phase-sensitive quantities. However, to incorporate the pairing field in our simulations
requires a number of important considerations. Since our modeling assumes contact interactions,numerical studies have to account for the potentially divergent nature of the pairing field at both
short and long length scales by appropriate renormalization of the scattering potential. We demonstrate how to renormalize the scattering potential when momentum is represented on a discrete
grid. Furthermore, when solving for the initial condition of the system, instead of using an approximation that ignores the pairing field [64] in order to remedy issues associated with the gaplessenergy spectrum, we take an alternative approach in which the condensate, depletion and pairing
field are accounted for and we solve for the HFB theory self-consistently.
We investigate the metrological use of the atom pairs generated from the periodic driving
protocol. We analyze the Hamiltonian for the excitation fields and associate it with that of opticalsqueezing. The analogy implies that the variance in of the the quadratures of the excitation field
can go beyond the standard quantum limit. Since at high momentum, the quasiparticle basis is37
close to the momentum basis, we point out that the atoms are ideal candidates for the momentum
superposition state used in matterwave interferometry.
The chapter is outlined as follows. We present the model in Section 4.2 and provide details of the renormalization process in Section 4.3. At first we limit our discussion to the moststraightforward case of quasi-1D systems. In Section 4.4, we outline the numerical procedures necessary to obtain a self-consistent ground state solution to the HFB theory for a weakly-interacting
trapped quantum gas, and quantify the quantum depletion as well as the pairing field amplitude.
In Sections 4.5 and 4.6, we use the time-dependent HFB theory to show that the modificationof the interaction strength through modulation of the scattering length parametrically amplifies acertain quasiparticle mode and generates a matter-wave solution that is analogous to a squeezed
state of light. In Section 4.7, we use these results to explore the possibility of future experimentsthat utilize interferometry to probe the pair correlation amplitude. We consider two methods thatcreate a phase difference between the driving field and the pairing field, and consequently lead tothe possibility for constructive and destructive interference in the matter-wave density. Finally, in
Section 4.8, we extend the results to quasi-2D so that they can be compared with the experimental
observations of angular correlations in the firework pattern, where the atoms were confined in apancake-shaped confining potential well and were ballistically expanded.

4.2

Formalism: General many-body field theory
We begin from the many-body Hamiltonian that describes a weakly interacting Bose gas withpairwise contact interactions:
ˆ
H =


ℏ2 2
d x ψ̂ (x) −∇ + Vext (x) ψ̂(x)2mˆV+d3 x ψ̂ † (x)ψ̂ † (x)ψ̂(x)ψ̂(x) ,
2
3

†



(4.1)

where m is the mass of the atom and Vext is the external trapping potential. The field operators,ψ̂(x) and ψ̂ † (x), are bosonic operators that annihilate and create particles and obey commutation
relations [ψ̂(x), ψ̂ † (x′ )] = δ(x − x′ ). The strength of the interaction potential, V , is related to the38s-wave scattering length, a, by V = T Γ , where T = 4πℏ2 a/m is the three-dimensional T -matrix(here it is actually a simple scalar and not a matrix since we consider the regime in which there isno dependence of the scattering phase shift on energy) and Γ is the dimensionless renormalization
factor that will be fully discussed in Section 4.3.
Since we intend to explore excitations from a BEC, we assume that the field operator iswell described by a mean field amplitude describing the atom condensate, ϕa (x), and a fluctuating
component, i.e.,ψ̂(x) = ⟨ψ̂(x)⟩ + δ ψ̂(x) = ϕa (x) + δ ψ̂(x) .

(4.2)

where ⟨δ ψ̂(x)⟩ = 0. The second-order terms—normal and anomalous densities—are defined respectively as,
GN (x, x′ ) = ⟨δ ψ̂ † (x′ )δ ψ̂(x)⟩ ,GA (x, x′ ) = ⟨δ ψ̂(x′ )δ ψ̂(x)⟩ .

(4.3)

Both of these play an important role in the dynamics of the non-condensate component of the system
we are interested in. In particular, the diagonal elements of the normal density, GN (x, x), represent
the physical non-condensate atom densities at position x and are therefore positive semi-definite.
The off-diagonal elements represent the matter-wave correlations of the non-condensate atoms thatare characterized by quantities such as the de Broglie wavelength and effective temperature. The
anomalous density, GA (x, x), is the pairing field that characterizes the two-particle correlations inthe system.
If we assume that the the field fluctuations are Gaussian, one can drop the third-order cumulants, and expand the fourth-order quantities in terms of the second-order cumulants when deriving
the evolution equations. In practice, this involves repeated application of Wick’s theorem [73]. The

39
resulting equations of motion are closed and can be written in detail for the condensate as;∂ϕa (x)
iℏ
∂t

=


ℏ2 2
−
∇ + Vext (x) ϕa (x)2m+ V |ϕa (x)|2 + 2GN (x, x) ϕa (x)+ V GA (x, x)ϕ∗a (x) ,(4.4)

for the normal density as;iℏ∂
GN (x, x′ ) = H′ (x)GN (x, x′ ) − H′ (x′ )GN (x, x′ )∂t+∆(x)G∗A (x, x′ ) − ∆∗ (x′ )GA (x, x′ )(4.5)

and for the anomalous density as;iℏ∂
GA (x, x′ )∂t

= H′ (x)GA (x, x′ ) + H′ (x′ )GA (x, x′ )+ ∆(x) G∗N (x, x′ ) + δ(x − x′ )+ ∆∗ (x′ )GN (x, x′ ) .

(4.6)

Here we have simplified the notation by introducing two energy functionals,
ℏ2 2∇ + Vext (x) + 2V |ϕa (x)|2 + GN (x, x) ,2m∆(x) = V ϕa (x)2 + GA (x, x) ,H′ (x) = −(4.7)

for the single-particle self-energy and the gap, respectively. Due to the fact that we neglect explicit
three-particle and higher correlations, the validity of this approach is restricted to the dilute gas
regime. Note that Eq. (4.4) can be simplified to the Gross-Pitaevskii equation if terms involving
GN (x, x) and GA (x, x) are dropped. In this case, the time-independent energy eigenvalue represents the chemical potential, µ, so that iℏϕ̇a = µϕa . The delta function in Eq. (4.6) arises from thebosonic commutation relation of the field operators and can therefore be interpreted as a quantum
effect. A number of quantities are conserved in this evolution; in particular, the total atom numberˆN =
dx |ϕa (x)|2 + GN (x, x)

(4.8)

40
is invariant under time evolution governed by Eqs. (4.4)–(4.7).
In order to see how the anomalous density is related to the vacuum pair wavefunction for theinteratomic separation of two atoms, we may neglect the mean-field density and the normal density
in Eq. (4.6), and then the eigenvalue equation is simplified to−ℏ2 2∇ GA (r) + V δ(r)GA (r) = 2µGA (r) ,
m

(4.9)

where r = x−x′ . Eq. (4.9) can be identified as a one-dimensional Schrödinger equation of a fictitious
particle of reduced mass m/2 scattering off a potential V δ(r). Then GA (r) is interpreted as theresulting eigenstate wavefunction corresponding to the familiar two-particle scattering solution of
the equation written in terms of the relative coordinate.

4.3

Renormalization of the Scattering Potential
The Dirac delta function in Eq. (4.9) implies that we are implicitly building a scatteringmodel from a contact interaction. This is convenient as it simplifies the resulting field theory,
but care must be taken to account for divergences that can arise at small and large scales. In
general, this is remedied by renormalization of the potential strength. In order to carry out thisrenormalization procedure, we begin from the formal scattering theory [91], where we define thebare scattering potential operator, V̂ , which has units of energy, and thereby expand the T -matrix
in an order-by-order series;
T̂= V̂ + V̂ G0 V̂ + V̂ G0 V̂ G0 V̂ + . . .
= V̂ + V̂ G0 T̂ .

(4.10)

Here G0 is the bare single particle propagator,
G0 =1
E − Ĥ0 + iϵ,

(4.11)

the scattering energy is E, the dispersion relation is Ĥ0 = p̂2 /(2m) with p̂ the momentum, and weneed to implicitly consider the limit ϵ → 0. The T -matrix elements are T = ⟨k′ | T̂ |k⟩, where |k⟩ is41
the wavenumber basis state. For the low energy scattering limit, the T -matrix becomes independentof E, and does not depend on k or k′ . In this case the T -matrix is well characterized by a constantscalar associated with the s-wave scattering length, as mentioned earlier, i.e., T = 4πℏ2 a/m.
Further considerations have to be made when one or more dimensions are effectively frozen
out due to imposing a strong confining potential in these dimensions. Without loss of generality,
let us consider the strong confining potential to be a harmonic potential with oscillator length
given by l⊥ . If one dimension is frozen out, an effective quasi-2D geometry is realized, and if two
dimensions are frozen out, an effective quasi-1D system is generated. If we denote the number of freedimensions by n ∈ {1, 2, 3}, the appropriate T -matrix expression, Tn , for the reduced dimensional√case can be related recursively by T3 = T and Tn−1 = Tn /( 2πl⊥ ) [8][75].
The process of renormalization connects the T -matrix, Tn , to the strength of the potential, Vn ,by expanding Eq. (4.10) in the momentum basis, and this connection depends on the dimensionalityof the system,

ˆ K+
Tn = Vn + Vn
K−

Tn
dn k
,
n
(2π) E − ℏ2 k2(4.12)

2m

where the critical element here is the introduction of K− and K+ as infrared and ultraviolet
momentum cutoffs, respectively. The cutoffs have to be chosen from an appropriate asymptotic
limit in order to accurately capture the dynamics of interest. The renormalization procedure canbe represented by the introduction of a parameter, Γn , defined by solving Eq. (4.12) for Vn . This
gives the solution,
Vn =
where

Tn
≡ Tn Γn ,1 − αn Tnˆ K+

αn = −

K−

dn k
1
.
n
(2π) E − ℏ2 k2(4.13)(4.14)

2m

In order to illustrate the behavior of αn , we consider the solution to scattering equation, Eq. (4.9),
where the stationary energy eigenvalue is E = 2µ and the mass is replaced by the reduced massof two particles, m → m/2. Solving this system of equations has a character that depends onthe dimensionality. In three dimensions, we may set K− = 0, and perform the integral to give42
α3 = mK+ /(2π 2 ℏ2 ) for a particle scattering at low energy, E → 0 [46]. In 2D, the integral scales
logarithmically and has both ultraviolet and infrared divergences. In 1D, there is an infrared
divergence so K− must be non-zero but we may set K+ to infinity.
We do not provide all the details here, since, in practice, these are formal considerations
that do not actually affect our numerical simulations. Indeed there are actually no divergences
introduced that require the introduction of momentum cutoffs to rectify when the momenta are restricted to values on discrete and finite grids. This is always the case in a numerical computer modelthat aims to describe a realistic experiment. In such a discrete representation of possible momenta,
it is preferable to simply calculate a finite sum over a specific partition instead of evaluating the
−1
P −1 2 2ℏ ki /m − 2µ ,continuous integral analytically. This implies a numerical evaluation of N
i=0
giving αn , and therefore determining Vn for a given T -matrix, which replaces V in the HFB equations, i.e., Eqs. (4.4)–(4.6). Here the subscripts i label individual discrete momenta, and thus {ki }represents the momentum grid, with N is the total number of grid points.
We carry out this renormalization procedure for all the results that we present in this chapter.
For each calculation, we verify that the numerical results are independent of the details of the
momentum grid on which the field theory is represented.

4.4

Self-consistent ground state solution
In order to find a self-consistent solution to prepare an initial condition for the subsequent

time evolution, the first step will be to consider the non-condensate component to be absent,
and to find a ground state representation of the condensate by solving the GPE. We then use
this condensate field as input into the time-independent equations for the normal and anomalous
densities, and diagonalize the resulting HFB self-energy matrix to find the quasiparticle basis. As
we will see, this process exhibits a defect in the zero-energy subspace (i.e., the eigenvectors do
not span the space). The intepretation is that the eigensolution is not stationary and cannot be
used as an accurate description of the initial condition for subsequent time evolution. We therefore
reintroduce the non-condensate terms that we have just found into the equations for the condensate,43
normal density, and anomalous density and solve again the system of equations, giving rise to an
iterative method that generates an accurate self-consistent initial condition.
Our approach will be to begin by first fully describing the necessary procedure using thesimple case of quasi-1D where the problem is most easily tractable. However, higher dimensions
can be treated in a similar method to the manner we present (we will consider quasi-2D laterin Section 4.8). The reduction to one-dimensional behavior requires the transverse confinement
conditiona2
n1D l⊥≪1(4.15)

to be satisfied, where n1D is the one-dimensional density [68][51], and as defined previously, l⊥ isthe harmonic oscillator length in the two strongly confining directions, here assumed to be equal.
The first part of our numerical algorithm is to solve for the ground state of the GPE,
∂ϕa (x)
iℏ
=
∂t




ℏ2 2
−
∇ + Vext (x) ϕa (x) + T1 |ϕa (x)|2 ϕa (x) .
2m(4.16)We use imaginary-time propagation to derive the lowest energy solution and its energy eigenvalueµ representing the associated chemical potential. Then, the mean field solution, ϕa (x), can be used
as a parameter to construct the self-energy matrix;


ΣA 
 ΣN
Σ=
,
∗ −Σ ∗−ΣA
N(4.17)

where
ΣN

= −

ℏ2 2
∇ − µ + 2T1 |ϕa (x)|2 ,2m

ΣA = T1 ϕa (x)2 .

(4.18)The self-energy matrix has dimensionality 2N × 2N where N is the size of the single particle
basis, as defined previously. This energy operator is most simply expressed in the position basis,
where x ∈ (0, L], since in that representation the potential terms including the mean-field appear
as diagonal blocks. The eigenstates of Σ are the Bogoliubov quasiparticles. Since the matrix
satisfies σz Σ = Σ † σz , where σz = diag(IN ×N , −IN ×N ), the eigenenergies come in pairs of positive

44
and negative values, ±ϵk , and the corresponding eigenstates are wk = (uk (x), vk (x))T and w−k =(vk∗ (x), u∗k (x))T . The eigenstates are normalized by satisfying the constraint
ˆ L
N
X


L
dx |uk (x)|2 − |vk (x)|2 −→|uk (jL/N )|2 − |vk (jL/N )|2 = 1 ,N
0(4.19)

j=1

Although this construction may appear standard and straightforward, there is a well-known andimplicit subtlety when examining the solutions to this eigensystem. When investigating the zeroenergy eigensolutions, one finds a pair of eigenstates that are colinear (equal up to a multiplicative√√
scalar) that have the form P = (ϕa (x)/ 2N , −ϕ∗a (x)/ 2N )T . This solution creates two significant
issues. First, the colinear eigenstates cannot be normalized by Eq. (4.19). Second, they do notspan the two-dimensional subspace of the Hilbert space corresponding to zero-energy.
The origin of this mathematical fact has an intuitive explanation. It arises from the approximations that lead to this self-energy matrix, that is, by fixing the condensate solution as an
unchanging parameter, one builds an unphysical model that implicitly allows the unconstrained
growth of a zero energy mode as a function of increasing time. Consequently there is no stationary
solution. This has to be remedied, for example, through a self-consistent approach in which thecondensate is treated as a variational parameter, in order to allow us to extend the formalism sothat it may be applied to our system of interest.
We begin by determining the remaining eigenvector to fully span the zero-energy subspaceby employing the Gram-Schmidt orthogonalization method to numerically calculate the remaining
basis vector. In this way we determine an eigenvector solution Q = (q(x), −q ∗ (x))T such thatL †L †N Q σz wk = 0 for all k ̸= 0, and normalize it to N Q σz P = i [97]. The addition of this vector to

the eigenvectors of the self-energy completes the basis of the vector space. The reason that this isimportant is that it allows the field operator to be expanded asδ ψ̂(x) =N−1X

ϕa (x)
θ̂ + iq(x)L̂(uk (x)b̂k + vk∗ (x)b̂†k ) − i √2N
k=1(4.20)

where b̂k and b̂†k are bosonic annihilation and creation operators for the quasiparticles. We haveintroduced θ̂ and L̂ as a pair of canonically conjugate operators that fully describe the zero-energy
mode and obey the cannonical commutation relation [θ̂, L̂] = i.

45
It is convenient to identify two special combinations of P and Q in order to give a concise√
expression for the completeness relation. We define w± = (∓P + iQ)/ 2, along with the matrix

W =



w+ , w1 , ..., wN −1 , w− , w−1 , ..., w−(N −1)(4.21)

so that the following completeness relation is satisfied;L †W σz W = σz .
N

(4.22)

This allows the particle annihilation operator to be written asδ ψ̂(x) =X(uk (x)b̂k + vk∗ (x)b̂†k ) ,(4.23)

k∈S

where the sum is over the elements of the index set S = {+, 1, 2, . . . , N − 1}, and b̂+ is the√annihilation operator for the zero-energy mode given by b̂+ = (iθ̂ + L̂)/ 2.
At this point, we have determined the quasiparticle basis, and can populate that basis witha given set of probabilities in order to generate particle distributions. In particular, we would liketo derive the normal GN (x, x′ ) and anomalous GA (x, x′ ) densities that are essential elements of the
HFB theory. To begin with we construct the Hermitian density matrix:′′GA (x, x )GN (x, x )
G = 

′
′
′
∗
∗
GA (x, x ) δ(x − x ) + GN (x, x )= W ΠW † ,(4.24)

where the population matrix Π has the form



p
q


Π=
.
q∗ I + p(4.25)The diagonal elements of p are the populations of each quasiparticle ⟨b̂†k b̂k ⟩, and the off-diagonal
elements represent the correlations between different quasiparticles. In the ground state, p = 0 and
q = 0. The identity, I, on the lower-right block is interpreted as a bosonic analog to the Dirac sea[33], in which the negative energy states are occupied by boson holes. When there is an excitation,

46
a pair of one particle and one hole is created, and therefore p appears in both the upper-left andthe lower-right block, as shown in Eq. (4.25).
This formalism now allows an extremely concise representation of the full dynamical evolution
encapsulated in Eqs. (4.5) and (4.6);iℏ∂G= ΣG − GΣ †∂t

(4.26)

where G is defined according to Eq. (4.24). The consequence of completing the basis by establishing
the missing eigenvector through Gram-Schmidt orthogonalization is now evident. If we begin withthe bare Σ, as defined in Eq. (4.17) and initialize G to the ground state (meaning p = 0 and q = 0)of the corresponding eigenbasis, then when Eq. (4.26) is propagated from this initial condition, it
is evident that the solution is not stationary. The number of non-condensate atoms is seen to growas ∼ t2 , as shown in Fig. 4.1. This implies that we have not in fact determined the correct groundstate.
This problem arises because, using the language of quantum optics, we are effectively assuming that the condensate is a coherent field that may act as an infinite classical pump and canprovide a reservoir source for introducing an infinite number of atom-pairs. Furthermore, it does
not cost any energy to introduce a zero-energy quasiparticle within this framework. This is clearly
unphysical for a number of reasons including the fact that, as can be seen in Eq. (4.4), the factorof two in front of the interaction between the condensate and the non-condensate atoms meansthat it actually costs energy to take away atoms from the condensate and move them into the noncondensate fraction, providing the interactions are repulsive (scattering length positive). There issome literature that suggests simply dropping the zero-modes entirely to remedy this problem, forexample, Ref. [9]; however, this generally violates the fundamental commutation relations of thebosonic field operator and therefore the uncertainty principle, so we do not employ that approachhere.
We instead employ an alternative solution by including the second-order terms to generalize
the self-energy matrix. This means that we modify Eq. (4.18) to include the effects of the normal

47


´
Figure 4.1: Quantum depletion N1 GN (x, x) dx as a function of time (the proportion of
non-condensate atoms at zero-temperature) simulated with the gapless HFB theory (i.e., using
Eqs.(4.16), (4.18), (4.21), (4.24), and (4.26)) gives a depletion proportion that initially scalesas ∼ t2 .

48
and anomalous densities, and then introduce the renormalization of the T -matrix to give
ΣN

= −

ℏ2 2
∇ − µ + 2V1 [|ϕa (x)|2 + GN (x, x)] ,2m

ΣA = V1 [ϕa (x)2 + GA (x, x)] .

(4.27)

with both Eq. (4.17) and Eq. (4.26) unmodified. In order to be consistent, however, we must also
generalize the GPE, Eq. (4.16), toiℏ∂ϕa (x)
∂t

= −

ℏ2 2
∇ ϕa (x) − µϕa (x) + V1 [|ϕa (x)|2 + 2GN (x, x)]ϕa (x)2m+ GA (x, x)ϕ∗a (x, x) .

(4.28)

Note that the ground state solution of the GPE is stationary, and thereby determines the valueof the chemical potential that enters the renormalization (see Section 4.3). Since GN (x, x) andGA (x, x) are functionally dependent on the eigenstates themselves, the problem is nonlinear, and
it is necessary to solve the generalized self-energy, Eq. (4.27), and the generalized GPE, Eq. (4.28),
iteratively until the equations are self-consistent [31]. We point out that this iterative process will
typically create a small gap in the energy spectrum of the system around zero energy, and the
problem of the unphysical non-stationary eigensolution that is caused by the zero-energy subspace
is no longer present. The resulting self-consistent solution is stationary under the evolution given
by Eq. (4.26) and provides an accurate ground state initial condition for the subsequent timedependent simulations that follow.

4.5

Dynamics of the time-dependent HFB system
In the experiment by Clark et al. [18], an external sinusoidally oscillating magnetic field isapplied, and therefore the scattering length is modulated in the forma(t) = adc + aac sin ωt,(4.29)

where adc is the initial scattering length, and aac is the amplitude of the oscillating component of
the scattering length at angular frequency ω. The dynamics of the system under this modulation

49

Figure 4.2: Solutions for a system with total atom number N = 6×105 in a 1D infinite potential well
2 /L. (a) Ground stateof size L, with the scattering potential between atoms given by a = 10−4 l⊥condensate density found from the gapped self-consistent generalized GPE theory (i.e., replacing
Eq. (4.16) with Eq. (4.28) and Eq. (4.18) with Eq. (4.27)). The length scale over which thecondensate density falls to zero at the edges of the box is known is the healing length. (b) Solution
to the normal density GN (x, x′ ) in the ground state as found from the self-consistent HFB theory.
(c) Solution to the absolute value of the anomalous density |GA (x, x′ )| in the ground state as foundfrom the self-consistent HFB theory.

50
is interesting to consider because the oscillating external field will inject energy into the system,and this will result in exciting atoms from the ground state into higher quasiparticle levels.
We begin our simulations by preparing the system in the self-consistent ground-state of the
HFB theory for a small positive value of the scattering length using the procedure just described. An
illustration of the resulting condensate, normal and anomalous densities are shown in figure 4.2.
After preparing the system in the ground state, we solve Eq. (4.17) and Eq. (4.26) using thegeneralized Eqs. (4.27) and (4.28) with a sinusoidal modulation of the scattering potential, i.e.,
V → V (t) = Vdc + Vac sin ωt .

(4.30)In order to interpret our results, we display the occupation probabilities via the projectionof G(t) onto the initial quasiparticle basis found from the self-consistent HFB Hamiltonian attime t = 0. The procedure is as follows. Since the quasiparticle eigenbasis matrix, W , satisfies
L†
N W σz W = σz , we may write

W −1 =L
σz W † σz .
N

(4.31)

Then, according to Eq. (4.24),
Π = W0−1 GW0†=

−1(4.32)L2
σz W0† σz Gσz W0 σz
N2where W0 is the original self-consistent quasiparticle basis determined for the initial condition. The
resulting population is shown in figure (4.3). The height of the peak in the off-diagonal block (i.e.,q) is notable since the coherence saturates the upper bound of the Cauchy-Schwartz inequality,δ(x − x′ ) + GN (x, x) GN (x′ , x′ ) ≥ |GA (x, x′ )|2 ,(4.33)

which in turn can be interpreted as confirming that the process of exciting quasiparticles fromthe condensate is maximally coherent. The diagonal elements, pk , can be measured by time offlight, since the quasiparticles transform into regular particles that can be detected during ballistic
expansion. In other words, when the kinetic energy greatly exceedes the interaction energy, the51

Figure 4.3: Snapshots of (a) absolute value of the matrix elements of the upper-left block of the
population matrix, |p|, and (b) absolute value of the matrix elements of the upper-right block of the
population matrix, |q|, at t = 0.1 (mL2 /ℏ), starting from the initial condition shown in figure 4.2 and2 /L and frequency ω = 1000 (ℏ/mL2 ). This
then continuously driven with amplitude aac = 10−4 l⊥frequency resonates with the quasiparticles with energy ϵk∗ = 500 (ℏ2 /mL2 ), corresponding to theresonant wavenumbers shown for reference as white lines (at k-index (πk)2 ≈ 500). At the resonantquasiparticle excitation a clear spike is evident. (c) Density profile of the condensate, revealingin general form the spatial dependence of the eigenmode function of the resonant quasiparticle
excitation.

52
k’s then effectively label the free momentum, i.e., kℏπ/L. As shown in figure (4.3), when the
periodic drive is turned on continuously for many cycles, essentially only one quasiparticle mode isresonantly amplified. That is consistent with the narrow spectrum. This physical process can beinterpreted as being due, as a consequence of the oscillating drive, to a photon with energy ℏω being
absorbed by a pair of atoms, with each of them getting half the energy, ϵk = ℏω/2. In addition, thephonon-like collective excitations that correspond to the observed wave-like patterns seen in thecondensate density can be interpreted as the Faraday patterns that typically manifest in different
kinds of parametrically driven fluids [88]. The pattern resembles the wavefunction density for the
single quasiparticle mode on resonance. This simulation illustrates that by careful engineering of
the drive, one can potentially prepare a variety of quantum states, selectively exciting atoms fromthe condensate field. We now show a few illustrative examples of interesting cases that employ thistechnique.

4.6

Dynamically Generating Squeezed Quasiparticle States
A squeezed state refers to a quantum state that has a reduced uncertainty in one degree of

freedom (‘squeezed’) at the expense of increased uncertainty in a canonically conjugate variable [44,104, 96, 105, 80, 81, 84]. Such states have been extensively studied in quantum optics and atomic
physics due to their utility in quantum metrology for producing measurement precision that exceeds
the limits derived from classical states. Here we will show how to use the resonant quasiparticle
excitation in order to generate a squeezed matter-wave state, anticipating that this could potentiallybe applied to quantum matter-wave interferometry.
By driving the system resonantly, we are effectively producing resonant pairs with well definedenergy, and this is reminiscent of nonlinear optical devices that down-convert pump photons into
signal and idler pairs. Here we will demonstrate that this correspondence is robust and quantitativeby demonstrating how one may calculate the squeezing parameter associated with the analogous
quantity that is regularly computed in the quantum description of light.
In order to do this we assume a weak excitation limit, so that GN (x, x) and GA (x, x) are small

53
compared to |ϕa (x)|2 . Furthermore, we consider the kinetic energy term in the time-dependent GPE
to be small, and then we can find a general solution for the condensate that has the formϕa (t) = ϕ0 e

iA cos(ωt)

= ϕ0

n=∞
X

Jn (A)einωt ,(4.34)n=−∞where A = Vac |ϕa |2 /ℏω and Jn (. . .) is the Bessel function of the first kind. We will limit our
discussion to the case of high modulation frequency, in which the photon energy associated withthe drive, ℏω, greatly exceeds the mean field shift associated with the drive amplitude, Vac |ϕa |2 , sothat A ≪ 1. In this case the n = 0 term completely dominates the series expansion and we candrop all other terms.
The initial stationary Hamiltonian for the fluctuations can be written as H0 =†k ϵk b̂k b̂k ,Pwhere ϵk is the energy of the k-th quasiparticle, and the transformation to a rotating frame involves
making the replacement of the quasiparticle operatorsb̂k → b̂k eiϵk t/ℏ .

(4.35)

The contact interaction term in the Hamiltonian can be derived from the interaction term ofEq. (4.1),
VacHI =sin ωt2ˆdx ψ̂ † (x)ψ̂ † (x)ψ̂(x)ψ̂(x)(4.36)

From this point, we keep only the second-order terms in δ ψ̂, because these terms correspond toexponential growth and therefore dominate the solution. In order to simplify the problem further,
we assume that the drive frequency ω corresponds to the resonance condition ω = 2ϵk /ℏ, andintroduce the rotating wave approximation, which allows us to keep only terms with e±i(ω−2ϵk /ℏ)t .

54
By representing δ ψ̂ in the quasiparticle basis, we obtain an effective interaction Hamiltonian
HI

≈
≈

ˆ


Vac
sin ωt dx 4|ϕa |2 δ ψ̂ † δ ψ̂ + ϕ2a δ ψ̂ † δ ψ̂ † + ϕ∗a 2 δ ψ̂δ ψ̂
2

ˆ
Vac
sin ωt dx 4|ϕ0 |2
2
X

X
×
u∗k′ b̂†k′ eiϵk′ t/ℏ + vk′ b̂k′ e−iϵk′ t/ℏ
uk′′ b̂k′′ e−iϵk′′ t/ℏ + vk∗′′ b̂†k′′ eiϵk′′ t/ℏ
k′

+ϕ20

k′′

X

u∗k′ b̂†k′ eiϵk′ t/ℏ + vk′ b̂k′ e−iϵk′ t/ℏ

X

k′

+ϕ∗0

X
2

u∗k′′ b̂†k′′ eiϵk′′ t/ℏ + vk′′ b̂k′′ e−iϵk′′ t/ℏk′′uk′ b̂k′ e−iϵk′ t/ℏk′

+ vk∗′ b̂†k′ eiϵk′ t/ℏX

−iϵk′′ t/ℏ

uk′′ b̂k′′ e+ vk∗′′ b̂†k′′ eiϵk′′ t/ℏ



k′′

ˆ


Vac
≈
dx 4|ϕ0 |2 vk uk + ϕ20 vk2 + ϕ∗0 2 u2k b̂k b̂k
4i

ˆ


− dx 4|ϕ0 |2 u∗k vk∗ + ϕ20 u∗k 2 + ϕ∗0 2 vk∗ 2 b̂†k b̂†k(4.37)

This corresponds to the interaction Hamiltonian of a parametric amplifier in nonlinear quantum optics, namely HI = −iℏ χ2 (â2 − â†2 ), where χ represents the second-order nonlinear susceptibility that corresponds to the squeezing rate. We refer to the resulting time-evolved state as asqueezed quasiparticle state since the analog is an archetypal system for creating squeezed statesof light. This mapping allows us to extract the squeezing rate, i.e.,
χ=

Vac
2ℏ

ˆ



dx 4|ϕa |2 vk uk + ϕ2a vk2 + ϕ∗a 2 u2k ,(4.38)and the squeezing parameter increases with time at this rate, i.e, ξ = χt. If we choose the phasesof ϕa , uk and vk appropriately, then χ is real. As expected from the known optical solutions, thepopulation in the k-th quasiparticle mode grows proportional to sinh2 (χt). Figure 4.4 shows the
population as a function of time at different modulation amplitudes. Since sinh2 (χt) → e2χt /4 atlarge t, one can extract the squeezing rate from the asymptotic slope of log pkk . We confirm thatthe squeezing rate is proportional to the modulation amplitude, as indicated by Eq. (4.38).
Squeezed states are characterized by reduced variance in one quadrature at the expense ofincreased variance in the other quadrature perpendicular to it. We define the quadrature for the
resonant quasiparticles asXθ ≡ b̂†k eiθ + b̂k e−iθ ,(4.39)

55

Figure 4.4: Population in the resonant quasiparticle mode (labelled k) driven with frequency ω =5500 (ℏ/mL2 ) as a function of time for different modulation amplitudes. The amplitudes are Vref(yellow), 2Vref (red), and 4Vref (blue), where Vref = 1.25 × 10−4 (ℏ2 /mL). The slope of log pkk atlarge time is equal to twice the squeezing rate, which is proportional to the modulation amplitude
as shown in Eq. (4.38). The squeezing rates calculated from the slopes of the curves in the intervalt = [0.3, 0.4] (mL2 /ℏ) are 4.3, 7.2, 14.1 (ℏ/mL2 ), and the squeezing rates calculated from Eq. (4.38)
are 3.6, 7.1, 14.3 (ℏ/mL2 ) respectively for amplitudes Vref , 2Vref , 4Vref .

56
where θ is the angle of the orientation of the quadrature. Then the variance is(∆Xθ )

2=b̂†k eiθ + b̂k e−iθ

2 

D
E2
− b̂†k eiθ + b̂k e−iθ(4.40)= 2pkk + 1 + 2Re{qkk e−i2θ } .
The variance as a function of θ is shown in figure 4.5. We see that the variance at certain quadraturephase angles, θ, of states produced by modulation of the scattering potential can fall below thestandard quantum limit. The standard quantum limit is the level generated by the uncertainty
principle under the assumption that the variance in all angles θ is uniform.
Although direct measurement of the squeezing may not be as straightforward to implement
as in its optics counterpart, it may be possible to observe directly the atom coincidence (since
the particles are produced in pairs) on detectors placed in directions corresponding to opposite
momenta, and thereby measure the second-order coherence. The direct analogue of phase sensitive
photodetection (homodyne and heterodyne detection, for example) is generally more complicatedto implement with atoms than light, but in the next section we propose a possible experimentthat could be used to perform an analogue of such interference measurements on the squeezed
quasiparticle distributions that are generated.

4.7

Interferometry with squeezed quasiparticles
In principle, the diagonal elements of the normal density are the quantities that can be directlyprobed with standard atomic density images, for example in dispersive, absorption, or fluorescenceimaging techniques. On the other hand, the off-diagonal elements of the normal density and theanomalous density cannot be directly observed since they are phase dependent quantities and havecomplex values that require an interferometric method to determine the phases. We investigate thephase dependence of the quasiparticle production by analysing two distinct methodologies. One
approach is a potential experiment that is capable of performing the phase measurement throughthe use of a protocol that is based on the Ramsey sequence widely used in atomic physics [45].
A second alternative approach is closely associated with a recent experiment by Hu et al. [36],57

Figure 4.5: The variance of the quadrature as a function of the angle, at t = 0.05 (mL2 /ℏ) withmodulation amplitude Vac = 2 × 10−3 (ℏ2 /mL), evaluated using Eq. (4.40). The dotted line is
the standard quantum limit, where the variance is equal to 1. For a certain range of angles, thevariance falls below the standard quantum limit. More specifically, at θ = 0.88π the variance has
minimum, which means measurements of the quadrature along this direction will have the greatest
precision.

58
who demonstrated that applying a phase shift to the oscillatory field after driving the system for aperiod of time will suppress the non-condensate atom number, and that a π phase shift results in
the greatest suppression.
Our Ramsey protocol is as follows. First, we apply a non-zero Vac for a period of time τto implement the first oscillatory field in the Ramsey sequence. We then set Vac zero for a brief
waiting period of time ∆t. During this interval the anomalous density evolves freely at the resonancefrequency, 2ϵk /ℏ, and because there is no external work done on the system, the number of noncondensate atoms remains essentially constant. Next, Vac is set to the same nonzero value as earlier
to implement the second oscillatory field, again for the same period of time τ . This sequence isillustated in figure 4.6(a). From our simulation results, shown in figure 4.6(b), we observe that the
number of non-condensate atoms oscillates as a function of the free evolution time, ∆t. This isbecause a phase difference θ = ω∆t accumulates between the anomalous density and the driving
field during the free evolution period. We account for this behavior by showing that the oscillations
observed are a consequence of the driving field in the second zone either amplifying or attenuating
the anomalous density depending on the accumulated relative phase.
For comparison, we now examine an abrupt phase change protocol based on the Hu et al.
[36] experiment. We consider the effect of the phase shift by first modulating the interaction fora period of time τ , then applying a phase shift θ to the oscillating drive, and repeating again theinteraction for a period of time τ , as shown in figure 4.6(c). The result of the final non-condensate
atom number as a function of the phase shift is shown in figure (4.6)(d). It is interesting tocompare this protocol and the resulting fringe pattern to that found from the first method. The
explanation is that the two methods both operate in a manner that is analogous to a Mach-Zehnderinterferometer, where interference fringes are seen in the recombination of light propagating along
two paths as the relative accumulated phase is varied. In the first protocol that we have presented,
the phase is accumulated in the anomalous density, whereas in the second method, a direct phaseshift is applied to the external field. We have observed that both methods result in an interference
pattern with high visibility fringes that allow direct access for the observer to probe the phase

59
behavior.
The two methods, the complete Ramsey sequence or the abrupt intermediate phase shiftchange, can be understood in a similar formalism. Both the phase shift change and the Ramsey
wait-time effectively generate a phase shift in the direction of squeezing. This manifests as a change
in the phase of the squeezing rate, i.e., χ → χeiθ , and is associated with the resonant quasiparticle
state evolving under the unitary operator,
χ iθ

2

†2

Uθ (t) = e− 2 e (b̂k −b̂k )t(4.41)

during the subsequent time evolution period. In the Heisenberg picture, the time-evolved operatorb̂k for the quasiparticle at index k at the end of the sequence is therefore given byb̂k (2τ ) = Uθ† (τ )U0† (τ )b̂k (0)U0 (τ )Uθ (τ )ED
= cos2 χτ sin θ sinh2 χτ (1 + cos θ)b̂†k (2τ )b̂k (2τ )+ sin2 χτ sin θ cosh2 χτ (1 + cos θ) .

(4.42)

Note that at the special point θ = π, Uπ (t) = U0 (−t), and the second period of modulation
simply reverses the effect of the first period of modulation, so that the final population is zero.
However, we can see that in the numerical simulation, the final number of non-condensate atoms at
a phase shift of π is non-zero. This is because the analytic result is derived using the rotating-wave
approximation, and in the full simulation, the populations of the off-resonance quasiparticles are
not fully reversed due to the influence of the other terms that were dropped. In this case, thesecond period of the modulation may further increase their populations even at the special point,θ = π, leading to the observed finite non-condensate population. As a consquence, the degree towhich the excitations can be fully reversed can be interpreted as a measurement of the fidelityof the protocol for producing quasiparticle squeezing. fidelity of the preparation of the squeezed
quasiparticle state.

60

Figure 4.6: The two methods for probing the phase of the quasiparticle squeezing. In both cases,
the scattering length at first oscillates at frequency ω = 5500 (ℏ/mL2 ) for a period of time τ . Att = τ , in (a) the amplitude of the oscillation remains at zero for a time interval ∆t, and then the
scattering length again oscillates for another period of time τ , and in (c) a phase shift θ is applied
to the oscillation . Panels (b) and (d) show the resulting non-condensate fraction at the final time
as a function of ∆t or phase shift θ, respectively, for the two cases.

61

4.8

Quasi-2D system
In order to make a more robust connection with the recent Bose firework experiment [18],we would like to generalize the formalism we have presented from a quasi-1D gas trapped in abox potential to a quasi-2D gas that is initially trapped by a circular potential with the third
out-of-plane direction frozen. Although this geometry adds new degrees of freedom to our previous
analysis, we may exploit the fact that the circular system possesses cylindrical symmetry, so that the
wavefunction of the condensate can be solved effectively as a 1D problem in the radial coordinate.
Note that the quasi-2D system differs from the quasi-1D system in a number of important ways. The
momentum correlations will manifest as angular correlations that may be detected by looking foratom-atom coincidence on two detectors aligned in opposite directions. Furthermore, the divergenceproperties of the renormalization problem are qualitatively different in two dimensions, as discussed
previously.
We begin by writing the fluctuations in the field operator in the quasiparticle basis using
appropriate indices for two dimensions,δ ψ̂(r, θ) =X∗(r)e−ilθ b̂†k,l ,uk,l (r)eilθ b̂k,l + vk,l(4.43)

k,lwhere k corresponds to the excitation in the radial coordinate and l represents the angular momentum quantum number. The angular momentum will modify the form of the kinetic energy
for the 2D quasiparticles by including a new centrifugal term, ℏ2 l2 /2mr2 , that arises physicallyfrom circulation about the trap center. Due to cylindrical symmetry, the normal and anomalous
densities should be functions of only three real variables, two radii and a relative angle, which wedenote by r1 , r2 , and ϕ ≡ θ2 − θ1 , respectively. For both normal and anomalous densities, we write
the functions in terms of their expansion in angular momentum,
GN,A (r1 , r2 , ϕ) =X

(l)

GN,A (r1 , r2 )eilϕ .

(4.44)l

The time evolution can then be solved by substituting this expansion into Eqs. (4.4)-(4.6) and using
the appropriate form for the two dimensional kinetic energy.

62
The first case we consider is for the situation in which the circular trap potential well isinfinite and has radius R0 ,Vext (r) =



 0

r < R0(4.45) ∞ otherwise
We prepare the quantum gas in the ground state with a small repulsive scattering length a in orderto stabilize the system mechanically. The repulsive interactions are characterized by the appropriate
2D T -matrix, as discussed in Section 4.3. Procedurally, we carry out a similar sequence of steps tothose previously discussed for quasi-1D. First we solve the GPE using imaginary-time propagation,and use that mean-field solution as the first iteration for the solution of the HFB equations, ignoring
the non-condensate terms in the HFB self-energy. As before, this solution is non-stationary andwe must iterate between the GPE and HFB solutions in order to find a self-consistent solution
whose resulting evolution gives rise to densities that do not depend on time. The resulting threecomponents, the condensate, the normal density, and the anomalous density, are illustrated infigure (4.7). Note that the anomalous density diverges in general as the Hankel function of thefirst kind as a function of the relative distance |r1 − r2 | close to the origin. This is an analytic
result that can be derived by solving the scattering equation Eq. (4.9) in 2D [3]. This emphasizes
an important point; the anomalous density cannot be accessed directly in experiment and does notform an observable.
Now that we have prepared an accurate initial state, we can then begin to examine its time
evolution when subjected to a drive via a modulation of the scattering length. Similar to what we
saw in quasi-1D, the modulation leads to excitation of quasiparticles with energies on resonancewith the modulation frequency. Figure (4.8) shows the normal and anomalous density as a functionof time and relative distance at the center of the trap. A principal feature of the radial density
dependence is the appearance of phonon-like excitations with well defined wave-number. The noncondensate density increases monotonically with time, as is consistent with the squeezing picture
discussed earlier. On the other hand the anomalous density oscillates in time tracking the external
field.

63

Figure 4.7: Ground state solutions for a system with total atom number N = 6 × 105 in an infinite
circular box of radius R0 . The scattering length is set to a = 3.99×10−5 l⊥ . (a) Condensate density
as a function of radial position. (b) Quantum depletion density as a function of radial position.
(c) Anomalous density with the center of mass position at the center of the trap i.e., r1 + r2 = 0,
as a function of the relative distance |r1 − r2 |. The divergence that scales as the Hankel functionof the first kind close to the origin is a result of 2D scattering theory.

64

Figure 4.8: Normal and anomalous densities in a system that is prepared in the ground state of
the self-consistent HFB solution with a positive scattering potential adc = 3.99 × 10−5 l⊥ , andthen subjected to the modulating drive with angular frequency ω = 1200 (ℏ/mR02 ) and constant
amplitude aac = 3.99 × 10−5 l⊥ . (a) Non-condensate density as a function of radial potitionr = |r1 + r2 |/2 (the origin is on the right) and time, i.e., GN (r, |r1 − r2 | = 0, t). (b) Magnitude of
the anomalous density as functions of the relative coordinate |r1 − r2 | (the origin is on the right)
and time, i.e., |GA (r = 0, |r1 − r2 |, t)|. Only a small time interval beginning at t = 0.6 (mR02 /ℏ) isshown, so that the oscillations are resolved.

65
In order to capture the the dynamics of the high momentum atoms emitted outwards in the
Bose fireworks experiment, we extend our simulations to a system with a finite trap potential thatis higher than the initial chemical potential but lower than the kinetic energy of the excited atoms,
using a smooth hypertangent functional of form,
Vwell
Vext (r) =2
r − R0 
1 + tanhζR0(4.46)

where Vwell and ζ are positive constants. The form of this external confining potential was chosen
to reduce numerical artifacts. Figure (4.9) shows snapshots of the condensate and non-condensate
densities at time t = 0.04 and 0.08 (mR02 /ℏ) during the drive. At t = 0.04 (mR02 /ℏ), we observethat the condensate density is pushed towards the edge of the trap, and that some non-condensate
atoms are generated. At t = 0.08 (mR02 /ℏ), we see phonon-like patterns in the condensate andnon-condensate densities that appear as ripples or waves. We see qualitatively a residual excitedcondensate that represents a component that does not have sufficient energy to overcome thepotential barrier, and a non-condensate density containing much more energetic atoms that isobserved to propagate outwards and leave the finite trap region. These are due to the energetic
quasiparticles created by the drive and form an experimentally observable quantity in ballistic
expansion images.
Figure (4.9) also shows the density currents, which indicate the flow of atoms, including thosein the condensate and the non-condensate. They are computed from
1d
J(r) = −
r dt

ˆ r
0

ρtot (r′ ) · r′ dr′ ,(4.47)

where ρtot is the total density. At t = 0.04 (mR02 /ℏ), the number of high momentum atoms is still
relatively small, and therefore the density current is also small. The density currents at r = 0.4 R0
point outwards while those at r = 0.8 R0 point inwards. Outside the trap there are few atoms, andthe currents are essentially zero. At t = 0.08 (mR02 /ℏ), the density currents near the edge of thetrap and outside the trap point outwards with large magnitudes, which represents the emission ofhigh momentum atom pairs.

66
Unlike in the experiment, where jet-like patterns were observed, our simulation results show
isotropic images. This is anticipated since the functionals we calculate for the condensate, normal,
and anomalous densities, represent probability densities and not individual realizations (i.e., they
are ensemble averages over many experimental realizations.) On the other hand, a given experiment is fundamentally different in that it represents a single trial that exhibits shot-to-shot noise
associated with the projection that occurs in a single quantum measurement. In order to model thisprojection noise it would be necessary to simulate quantum trajectories [34], rather than solving forthe density matrix evolution, and this may be done by adding white noise to the initial condensate
wavefunction (see for example Ref. [25]).
In order to demonstrate a quantitative comparison of the energy of the generated quasiparticles with respect to their ballistic motion, we present a numerical ‘time-of-flight’ calculation. In
this simulation, we measure the momentum of the atoms in the system by evaluating the speed
at which the gas expands. We define the effective size of the gas as the radius encircling a large
fixed fraction of the non-condensate (say more than 90%), so that at t = 0 we begin with size R0 .
Using this metric, from our simulations, we observe that initially a large amount of non-condensate
density is generated close to the center of the trap and most of the non-condensate atoms havenot left the finite trapping region, so that the size of the gas appears to be shrinking. However,later in the evolution and after a significant fraction of non-condensate atoms escapes the trap, theexpansion of the size of the gas becomes essentially ballistic (expansion size increasing linearly inptime), and the speed of the expansion is approximately (ℏω − 2µ)/m. The reason is as follows.
From the HFB equations, we know that



ℏ2 2
2
∇ − µ + 2V |ϕa | uk + V |ϕa |2 vk = ϵk uk .
−2m

(4.48)

Using the approximations V |ϕa |2 ≈ µ and vk ≈ 0 for large k, we may write−ℏ2 2∇ uk ≈ (ϵk − µ) uk .
2m

(4.49)

67
Futhermore, since the resonant quasiparticle energy is ℏω/2, we find1
mv 2 ≈
2



ℏω
−µ
2


.

(4.50)The speed of the expansion is evaluated from the slope of the fitted curve in figure 4.10. The fit only
includes data points from time t = 0.065 ∼ 0.08 (mR02 /ℏ) so that we avoid the initial transientswhere interaction energy within the condensate and between the condensate and noncondensate issignificant, and the motion is not ballistic. The speed we get from the slope agrees well with theanalytical result derived above.

4.9

Conclusion
In this chapter, we have developed a description of a condensate and non-condensate systemstarting from the many-body field theory Hamiltonian and deriving the evolution equations for thecondensate, normal density and anomalous density. Since we assumed contact interactions, thecontact potential may lead to divergences in the field theory at small and large momenta. We tookcare of this issue by properly renormalizing the scattering potential.
We solved the quantum fluctuations in the initial stationary state in 1D using the selfconsistent HFB theory, which does not involve any free parameters. Then, we simulated theamplification of the quantum fluctuations with a well-defined energy using the time-dependentequations. The amplification has aspects similar to the generation of squeezed states of light,
and we were able to verify that the variance of the quadrature phases can fall below the standard
quantum limit. We proposed to observe phase sensitive quantities through two alternate approaches
including Ramsey interferometry and discrete phase jumps. We showed how this is able to provide
information on the characterization of quasiparticle squeezed states. Finally, we showed simulation
results in 2D, and found that the excited non-condensate atoms eventually leave the trap andpropagate outwards at a well-defined speed, consistent with the experimentally observed time-offlight results. Although we showed simulation results for only quasi-1D and quasi-2D systems, 3D
systems would be interesting and can be analyzed systematically using similar approaches.

68

Figure 4.9: Densities of the condensate and the non-condensate at t = 0.04 and 0.08 (mR02 /ℏ).
The dark blue area is the range of simulation, and the white line indicates the trapping potential.
The yellow arrows represent the density currents, with their lengths proportional to the magnitude.
The system starts with the ground state with adc = 1.99 × 10−5 l⊥ and then the scattering length
is modulated with amplitude aac = 1.99 × 10−4 l⊥ for all time. The smoothing parameter of thefinite circular well was set to ζ = 0.2. At t = 0.04 (mR02 /ℏ), the atoms in the condensate arepushed towards the edge of the trap, but because they are of low energy, they cannot escape thetrap. The number of high momentum atoms is still relatively small, and therefore the density
current is also small. The density currents at r = 0.4 R0 point outwards while those at r = 0.8 R0
point inwards. Outside the trap, there are few atoms, and the currents are essentially zero. At
t = 0.08 (mR02 /ℏ), a large fraction of non-condensate atoms with high energy escape the trap,
and the density currents near the edge of the trap and outside the trap point outwards with large
magnitudes. For clarity, the density currents at t = 0.04 (mR02 /ℏ) are scaled up 3 times comparedto those at t = 0.08 (mR02 /ℏ).

69

Figure 4.10: The size of the non-condensate as a function of time. The size is defined as theradius R that encircles 96.2% of the non-condensate atoms. The system starts with initial chemical
potential µ = 26 (ℏ2 /mR02 ) and is driven with modulation frequency ω = 1000 (ℏ/mR02 ). The
speed of the propagation v is given
p by the slope of the fitted line (red), v = 30.8 (ℏ/mR0 ). Thisvalue is approximately equal to (ℏω − 2µ)/m, the speed of a particle whose kinetic energy is half
the photon energy minus the chemical potential. We have excluded the data points at times whenthe motion is not ballistic, including at early times where few non-condensate atoms escape thetrap and the interaction between the condensate and the non-condensate is significant.

70
Note that there is one other important consideration; our model does not include the collisions
terms in the kinetic theory that result in equilibration of the gas to its thermal state [53]. Neglecting
collisions is a good approximation for a dilute gas at low temperature, but implicitly requires us tolimit our discussion to the regime in which the time-scale between consecutive two-body collision
events greatly exceeds the time-scale of the quantum dynamics that we investigate.
We have demonstrated a method to generate momentum squeezed states that may be usefulfor metrology applications. This motivates us to further consider engineering the scattering length
as a function of time to generate two-mode squeezed states in quasimomentum that could be injected
into matter-wave interferometry. The entanglement properties of such states would be interesting
to investigate along with the metrological gain that arises from the quantum advantage.
The importance of pairing in this work also motivates us to consider a similar experiment on
fermions, where the interactions could be modulated by variation of the scattering length in the
BEC-BCS crossover regime. The motivation for this is simply that the pairing physics is closely
connected with the previously observed fermionic condensation. These considerations will be thesubject of future studies.

Chapter 5

Reinforcement-learning-based matter-wave interferometer in a shaken opticallattice

5.1

Introduction
Extensive efforts have been made to develop metrological devices based on interference anddetection of either electromagnetic waves or matter waves. By exploiting the quantum aspectsof superposition and entanglement, one can potentially achieve high sensitivity to phase shifts,
and this has inspired a wide variety of applications, including detecting gravitational waves [1],measuring the fine-structure constant [70], testing the universality of free fall [24], and inertial
sensing for GPS navigation [32, 76].
The usual direct-design approaches to engineering complex quantum systems that consistof many degrees of freedom or many particles are typically founded on experience with simpler
analogs or intuition for underlying mechanisms. This naturally leads to a paradigm that is mostaccessible in terms of understanding, but incorporates human bias that may potentially generate
non-optimal solutions. With this perspective in mind, we point out there are a few purely systematic
methods that are often used as a way to develop unbiased strategies, including optimal control [72]and optimization algorithms such as the Nelder-Mead simplex [63] or simulated annealing [43].
Utilizing these methods can allow one to explore solutions with more complicated forms with thepotential to reach closer-to-optimum control protocols.
Recently, it has been shown that quantum design may benefit tremendously from a branchof machine learning based on trial and error, known as reinforcement learning, that aims to employ

72
!$#

!"#
!"""#

!""#
!

!%#

%
$

!""#
!

!"#
!"#

!"""#
!"#

!"#

!""#

!"#

!""#

&

&

!"""#

Figure 5.1: Interferometers are composed of (i) a beam splitter, (ii) mirrors, and (iii) a recombiner.
Examples shown are; (a) an optical Mach-Zehnder interferometer (using half-silvered or conventional mirrors), (b) a Bragg interferometer (using three short light pulses of varying pulse area, π/2or π), and (c) a shaken lattice interferometer. In (c), the shaken lattice mimics the interferometercomponents by splitting, reflecting, and recombining the atoms, through design of a specific shaking
function for each case.

machines to find the optimal strategy for accomplishing a specific task. One of the reasons for its
success is that the learning framework is decoupled from human intuition, and therefore may explore
novel solutions that have been previously undiscovered. Moreover, in situations where problems
are so complex or extensive that naı̈ve brute-force algorithms are considered unfeasible, such as
playing games like go and chess [85, 87, 86], sophisticated reinforcement-learning algorithms havebeen developed to enable the machines to perform at a level that exceeds human capability. There
are many systems in quantum physics that fit naturally within this scope due to their underlying
complexity. For example, reinforcement learning has been applied to control and study phasetransitions of many-body quantum systems [14], to find strategies for quantum error correction[23, 62], to design quantum circuits [65], to prepare novel quantum states [13, 55, 58], to findprotocols for quantum communication [95], and to improve quantum sensors [83].
In this chapter, we will investigate the potential for reinforcement learning to improve on existing matter-wave interferometery. In Sec. 5.2, we describe the physical model for the interferometricsystem. We outline our reinforcement learning approach in Sec. 5.3 and explain the key concepts
crucial to understanding its application. In Sec. 5.4, we demonstrate how the reinforcement-learning
framework may be used to design a beam splitter and a mirror, and present the protocols that are73
learned. Following this we show how the components can be cascaded together to form a complete
interferometer. In Sec. 5.5, we evaluate the resulting performance of the interferometer through a
Bayesian statistical analysis.

5.2

Physical model
Although the topology of interferometers can vary somewhat, the archetypal design can bethought of as the Mach-Zehnder interferometer [54, 107]. This device is sensitive to the differential
phase accumulated between two alternate paths [see Fig. 5.1(a)]. Mach-Zehnder interferometers
are composed from three essential components: A beam splitter that separates the wave coherently
into two directions, mirrors that reflect both parts, and a recombiner where the waves are broughtback together and constructively or destructively interfere. While in the optical case, componentsthat split or reflect light beams are readily available (e.g., half-silvered or conventional mirrors), formatter-wave optics, the components have to be generated through the careful control of laser-atominteractions. The traditional way to do this is through Bragg diffraction [see Fig. 5.1(b)] [11, 40], inwhich a sequence of three short light pulses are used to separate, reflect and recombine the matterwave.
An alternative approach is possible for atoms moving in an optical lattice potential where
the intensity pattern can be shifted backward and forward in time [102, 100, 101]. The elementarycomponents, i.e., a beam splitter, reflectors, and a recombiner, in that case may be implemented
by “shaking” the lattice in a tailored pattern. If we denote the canonical position and momentum
operators of an atom by x̂ and p̂, respectively, with m the atom mass, the optical lattice system
can be described by the HamiltonianĤ(t) =



p̂2
V0
−
cos 2kL x̂ + ϕ(t) ,
2m
2

(5.1)

where kL is the laser wavenumber, and ϕ(t) is the time-dependent phase difference between two
counter-propagating lasers. The lattice, with constant amplitude V0 , is shaken through the variationof ϕ(t), since that parameter determines the position of the nodes and the anti-nodes of the standing

74
wave intensity.
Unlike a typical Mach-Zehnder interferometer, which is formally a two-port device and transforms quantum states according to simple SU(2) group rotations, the relevant eigenstates of the
shaken lattice potential consist of many Bloch states that can be coupled by the time-dependenceof the laser phase, ϕ(t). In some sense this situation represents a highly multi-path form of interferometry since the Bloch basis provides many accessible paths for the quantum wave function to
explore. While this establishes a rich evolution and could provide metrological benefit, it makes it
more difficult to design an intuitive control protocol. It is for this goal of obtaining a high performance solution among a complex landscape that we are led to consider reinforcement learning asa design approach.

5.3

Model-free reinforcement learning for design
In order to understand our design philosophy, we first need to address the important con-siderations to make when using machines to find solutions to design tasks. In this section we will
present the principal ideas and concepts that will form the underlying foundations of our learningbased methodology. Since there are a variety of methods available, we begin with a discussion of
the main structure that we will employ.
Reinforcement learning consists of a closed loop in which an agent invokes actions based on
the observed state of an environment, and an environment provides rewards to the agent based on
the observed outcome of its actions. The agent is tasked with the goal to discover the sequence ofactions for which it receives the highest possible terminal reward [90] (see Fig. 5.2). It does this bytrial and error, iteratively improving its actions in such a way as to corral the environment towarda target configuration.
We are primarily interested in a specific kind of reinforcement learning, known as modelfree learning [90], where the agent has no detailed insight as to the structure of the environment,
and cannot know a priori what effect an action will have on the environment state. Here theagent learns and makes decisions even when the environmental model may not be fully known.

75This general setting can even be applied to the situation in which the environment is represented
by an experimental apparatus that cannot be fully understood, and the reward is derived fromexperimental observation and feedback. Using this approach significantly expands the potential
scope of our work, since it is precisely this aspect that distinguishes model-free reinforcement
learning from classical optimal control methods. In classical optimal control, the optimization
depends heavily on the mathematical form of the model and, for complex systems, can be difficult
to implement in practice.
Typically, reinforcement learning problems are Markov decision processes, where the probability of transition to the next state is only dependent on the current state and the current action,and not on prior history. Evolution in this framework is referred to as a trajectory, where the stateis initially prepared, and then a sequence of actions and corresponding state updates are performed.
The trajectory steps continue until the chain comes to an end when a predetermined condition isreached, which might be a certain number of steps, or a terminating state of the environment. An
important concept is the idea of an optimum trajectory starting from any state, which is a trajectory such that the terminal reward is the maximum possible. We use what is known as Q-learning
to make the action decisions, as we now describe. Although this approach is not the only possible
choice, it is appropriate when the effect of an action on the environment state is deterministic,
which will be the case here.
Since there is potentially an enormous dimensionality of the state space and the action space,a brute-force search is not possible, and we employ a neural network to learn and to approximately
optimize the crucial decision-making task. That is, we define the input nodes of a neural network
to be a representation of an arbitrary vector in the state space. We ascribe to the output of theneural network a vector of quality factors, simplified to Q-values, Q(s, a) ∈ R, where each element
represents the desirability for a possible action, a, given an input state, s. Ideally the neural network
should take any arbitrary state as input and output a distribution in which the maximum Q-value
corresponds to the best next action to perform.
We train the neural network in the following manner. The agent will utilize a policy in order76!,%(*$!"#$% '())*" &$%#$%

φ

!"#$%

F

iψ̇t =
+ V̂ (φ)ψ
i∂|ψK̂ψ
= Ĥ(φ)|ψ&$'()*$+#$%-#./)0123#)'/%(*$

Figure 5.2: Framework of reinforcement learning. The agent chooses an action, the environment
responds to the action, and gives the next state and a reward as feedback. In our design task, theaction is the translation of the optical lattice, the environment evolves according to the Schrödinger
equation, the observation is the momentum population distribution, and the reward is a functionof the quantum fidelity. Illustrated on the left is an example neural network as the decision-making
agent. The neural network takes the state represented on its input layer, passes it through a hidden
layer, and generates a vector of Q-values at its output layer.

77
to decide on the next action for a given state, for example, by taking the action associated with the
maximum Q-value from the neural network output vector. The environment, in some initial states, receives the action a from the agent, and then updates the state, symbolically denoted asas−→ s′ .
The environment then reports back to the agent the reward, r(s′ ). As a step in a sequence, thisleads to the important concept of the discounted cumulative reward, known simply as the return.
The return, Y (s′ ), is defined as the combination of a current reward, r(s′ ), and every future rewardthat would be generated in an optimum trajectory seeded by s′ . Since any current action will havedecreasing influence on future decision making as the steps become more distant, it is useful todeweight each consecutive reward by a discount factor, γ ∈ [0, 1], i.e.,Y (s′ ) ≡ r(s′ ) + γr(s′′ ) + γ 2 r(s′′′ ) + . . .

(5.2)

where
a′

a′′

a′′′

s′ −
→ s′′ −→ s′′′ −−→ . . .
are steps along an optimum trajectory.
The learning process is now framed as the minimization of the squared difference between
the return, Y (s′ ), and the network output, Q(s, a), where the minimization is accomplished byvarying the internal weights and biases that constitute the neural network. Should the learning be
perfected, the squared difference would be zero, implying that the Q-value and corresponding returnwould be identical. We emphasize that, while mathematically complete, this learning definition isa formal construction that is intractable as written, since if we knew the optimum trajectory there
would be no need to carry out the learning at all. In fact, the only information that is accessible
about the future component of Y (s′ ) is the inference one can make from the neural network output
itself, introducing a self-consistency element to the optimization problem. In other words, the series
in Eq. (5.2) is approximated as
r(s′′ ) + γr(s′′′ ) + . . . ≈ max
Q(s′ , a′ )
′
a

(5.3)

78
so that
Y (s′ ) = r(s′ ) + γ max
Q(s′ , a′ ) .
′
a

(5.4)

The final result given in Eq. (5.4) encapsulates all the principal concepts that establish a closed
learning formulation that can be computationally implemented.
We employ several technical but important improvements that increase the convergence, suchas to replace the neural network used for the calculation of the future return, as approximated in
Eq. (5.3), with a second identical neural network, the target network. The target network is updated
at each step from the Q-network by a weighted average leading to increased stability [61, 94].
Another useful trick is to employ an ϵ-greedy policy in the training stage [90], which simply meansthat the action associated with the highest Q-network output value is not always chosen, but with ϵlikelihood, a random action is selected. This allows the neural network to explore a more extensive
state space and to avoid becoming trapped in a locally optimal solution. Finally, as is common
to many forms of machine learning, a variety of batch sampling methods, in this case for a replay
buffer of previously experienced transitions, (s, a, s′ ), can be used to assist in the efficiency andstability of the learning minimization stage [52]. For further details of the algorithmic structure we
refer the reader to Chapter 3.
The Q-network we will use for includes an input layer, one hidden layer, and an output layer.
Each layer is connected to the next one by a linear function followed by a ReLU (Rectified Linear
Unit) activation function, which represents the function max(x, 0) with the argument x being areal number. The relations between layers are explicitly written out as follows,
h = ReLU(W1 x + b1 )(5.5)

y = ReLU(W2 h + b2 ) ,(5.6)

where h denotes the hidden nodes, x denotes the input nodes, y denotes the output nodes, W1 , W2
are the matrices of weights, and b1 , b2 are vectors of biases. An example topology of the Q-network
is shown in Fig. 5.3.

79

!"#$%

'())*"&$%#$%

Figure 5.3: An example topology of the network that we used in the deep Q-algorithm.

80

5.4

Component design
We now use the Q-learning framework to design the essential elements for building an in-

terferometer, i.e. a splitter, reflectors, and a recombiner. One aspect that is worth emphasizingis that the standard Mach-Zehnder interferometer topology may not always be the optimal layoutfor determining an accumulated quantum phase. For example, it has been demonstrated in [57]that novel optical experimental setups that could potentially provide sensitivity surpassing that of
the Mach-Zehnder configuration can be designed with reinforcement learning. However, with theaim of simplicity, and even though it constrains our solution space, we will limit the discussion
here to the design of only the specified individual interferometer components. This will allow us toevaluate the efficacy of the learned protocols by direct comparison with the device characteristicsof the analogous Bragg interferometer. In the future, it may prove beneficial to allow for more
freedom in the design topology and to thereby ascertain whether this leads to further performancegains.
We have developed the learning terminology around words such as “state” and “environment”
where they have been given a specific meaning. We are now going to apply the methodology to thedesign task of controlling a complex quantum system. Unfortunately we will encounter the issue
that terms such as state and environment have conflicting and completely different meanings in thetheory of open quantum systems. In particular, to enable an efficient machine learning algorithmwith an addressable number of input nodes, the state space that we use is typically a compact
representation of a quantum state or operator. The environment, which denotes the system that
we wish to control and that is evolving under our imposed actions, is a completely different use of
the term to that common in open quantum system literature. In the next design descriptions, wewill provide the dictionary for mapping the learning terminology to the quantum variables.

81
5.4.1

Beam splitter

The beam splitter is a target-state design problem that can be specified as follows. We
initialize a quantum state in the ground state of an optical lattice potential, where the lattice
depth V0 is set to 10 Er (Er = ℏkL2 /2m is the recoil energy). The task is to shake the lattice insuch a way as to split the atoms from the ground state into an approximate equal superposition of|±p0 ⟩ momentum states. We will focus on p0 = 4ℏkL here, which we find to balance well the tradeoff between large momentum splitting and high-frequency components necessary for the shaking
function. If the target quantum state is a lattice eigenstate it is a stationary state under freeevolution, avoiding the need for a free evolution protocol. As a consequence, we assign the target
quantum state to be the third-excited Bloch state with zero quasimomentum. This eigenstate
has the attractive property that it closely approximates the desired superposition of momenta,
although not perfectly, and can be strongly coupled to the ground state by a reasonably simple
shaking solution.
The state space for learning, i.e., the vector space containing s, is a compact representation of
the Hilbert space. For this learning task, we define s to be a vector of populations (i.e., probability
distribution) for momenta chosen from a comb of discrete values, {2nℏk} for n ∈ Z, |n| < nmax .
This choice is motivated by the fact that the stimulated absorption and emission of photons only
couple momenta separated by this quantized spacing. The possible actions, a, are also reduced toa discrete set of possibilities, each action corresponding to a specific constant phase to be applied
during the corresponding step. This means that with these associations, the deep Q-network is
specified as mapping momentum populations at its input layer to Q-values at its output layer
associated with a finite set of possible lattice phases to be applied during the next update.
In order to apply the action generated by a choice of lattice phase and to thereby determine
the updated momentum distribution that results, we need to specify the environmental model. In
principle, this could be experimental or theoretical, and even if theoretical, could include stochastic
effects that arise from dissipation and noise. However, for simplicity, we will take the most basic

82
formulation. This means that our environmental model will update a pure quantum state throughthe Schrödinger equation evolution for an isolated system. The reward will be based on the fidelity
between the evolved quantum state and the target quantum state, i.e., the third-excited Blochstate, as calculated from the modulus-square inner product.
Actually, this environmental model allows us to accomplish two design tasks at the sametime. Once we find a protocol for splitting, the recombining task can be achieved by simply implementing the time-reversed protocol, a consequence of the time-reversal symmetry underpinning
the Schrödinger evolution and the time-reversal symmetry of both the initial and the target states.
A reinforcement learning outcome for the beam splitter is shown in Fig. 5.4, which illustrates
the learned shaking function, ϕ(t), and the corresponding evolution for the momentum probability
distribution that results. Even though this sequence of phases does not appear to have a predictablestructure, it is apparent that at the terminal time the momentum distribution has the anticipatedform of two well-defined peaks at ±4ℏkL . The calculated fidelity of the state is approximately
95%. It would be possible to optimize further by expanding the set of possible phase elements,
by increasing the total evolution time, or by optimizing the learning cycles. However, any realexperiment may possess imperfections and other aspects that are not well described by our isolated
system’s Schrödinger evolution, and the level of performance of our design is sufficient for the task
at hand.

5.4.2

Details of reinforcement learning for target state

Here we present how we formulate the control for splitting as a reinforcement learning problem. First of all, we define the input state to be a vector of the populations in the momentum
eigenstates,{|2nℏkL ⟩ | n = −3, −2, −1, 0, 1, 2, 3} .

(5.7)

83
The momentum states that are neglected here are barely populated the whole time. The actions
are chosen from a set of discrete phase values,{nπ/4 | n = −4, −2, −1, 0, 2} .
The reward is a function of fidelity, defined as




0
if d = 0
r =
,


F

 1−F
if d = 1
F= | ⟨ψtarget |ψ(t)⟩ |2 ,

(5.8)

(5.9)

(5.10)

where d = 1 when the maximum terminal time is reached or when the fidelity is higher than thethreshold 0.95, and d = 0 otherwise. The Q-network we use for this task has an input dimension of
7 and an output dimension of 6, as a result of (5.7) and (5.8). The neural network is implemented
in PyTorch. The values for the hyperparameters used for training this Q-network are listed inTable 5.1. The hyperparameters are determined from a randomized grid search carried out in theHyperparameters

Valuesγ
ταepisodesϵ decay
Hidden size
Batch size
Optimizer

0.999
0.999
0.001
20,0000.0001
98
64
Adam

Table 5.1: Choice of hyperparameters for the splitting task. The exploration probability ϵ decays
linearly for each episode at the rate represented by ϵ decay. The optimizer we use, ‘Adam’ (which
stands for adaptive momentum estimation), is a method for efficient stochastic gradient descent [42].

package Ray Tune, where a set of hyperparameters that results in the best reward at the end of
the training among the 20 sets is picked out.

84
5.4.3

Mirror

The second task is to reflect the momentum from |±p0 ⟩ to |∓p0 ⟩, that is, corresponding
to a matter-wave mirror. We should point out at the start that this task is essentially different
in character from the beam splitter, because it is not a target-state but a target-operator design
problem. The mirror is defined by the desired map,
α |p0 ⟩ + β |−p0 ⟩ → α |−p0 ⟩ + β |p0 ⟩ ,for any arbitrary α and β. This implies that the target unitary is any operator Ûtarget that satisfies

P̂ Ûtarget P̂ = |p0 ⟩ ⟨−p0 | + |−p0 ⟩ ⟨p0 | ,(5.11)

where P̂ = |p0 ⟩ ⟨p0 | + |−p0 ⟩ ⟨−p0 | is the projector onto the 2-dimensional subspace. The design
goal is, therefore, to shake the lattice as a function of time, such that the corresponding unitaryevolution, Û (t), approximates any Ûtarget that satisfies Eq. (5.11) at the terminal time. The vector
space containing s is defined as a compact representation, in this case, of the set of quantum
unitary operators. Given the mirror behavior, a suitable choice for s for this learning task is avector of norms of selected matrix elements of the unitary represented in momentum space, i.e.,| ⟨±p0 | Û (t) |p = 2nℏk⟩ |2 for n ∈ Z, |n| < nmax . For the set of possible actions, a, we use a different
parametrization of the phase than for the beam splitter design. We found that the quality of themirror is very sensitive to the frequency of the shaking function. The characteristic frequency,
12ωr , where ωr = ℏkL2 /2m is known as the recoil frequency of the lattice, corresponds to the kinetic
energy difference between |4ℏkL ⟩ and |2ℏkL ⟩, and works especially well. Consequently, we enforce
a sinusoidal ϕ(t) at frequency 12ωr , and at each half cycle, the agent chooses the action a as theamplitude of the sinusoidal shaking function selected from a small set of possible values. Thismeans here that the deep Q-network is specified as mapping a compact representation formed fromthe norm of a subset of elements of the unitary matrix to Q-values associated with the amplitudeof the sinusoidal shaking function for its next half cycle.

85
We also need to define what is meant by the environment in this case. The environmental
model establishes the unitary operator update during a shaking function half cycle. If, as before,
we take the environmental model to be a closed quantum system, the evolution is given by solving
the Schrödinger equation of motion for the time evolution operatoriℏdÛ (t)= Ĥ(t)Û (t) .
dt(5.12)The reward is associated with the channel fidelity in the ds = 2 dimensional subspace spanned by|±p0 ⟩ (see [71])F=1
ds (ds + 1)

h



i
Tr M M † + | Tr(M )|2(5.13)†
where the density operator is M = P Utarget
U (T )P .

An example outcome from reinforcement learning for the mirror is shown in Fig. 5.5. Thisillustrates the learned shaking function and corresponding momentum probability density that thisgenerates starting from −4ℏkL . From matrix elements of the resulting unitary operator, we also
verify that the phase relation between the coefficients of the ±4ℏkL states is preserved. As for the
beam splitter case, it would be difficult to anticipate the form of ϕ(t), and yet the resulting quantum
state evolution closely approximates that expected for a mirror. The channel fidelity reaches 93%
at the terminal time, which as discussed earlier could still be improved on, but nevertheless issufficient for our purpose.

5.4.4

Details of reinforcement learning for target unitary

The reflection operation is defined by a target unitary operator that satisfies Eq. (5.11). Our
goal is to control the phase as a function of time, such that the unitary operator

ˆ t

′
′
Û (t) = T exp−iĤ(ϕ(t ))dt /ℏ,(5.14)

0reaches the target operator at the terminal time. Here T represents the time-ordering operator
and Ĥ(ϕ) is defined in Eq. (5.1). To formulate the control problem as a reinforcement learning
problem, we assign the state to be the matrix elements of the unitary operator. Since before the86wave function enters the reflection stage, most of its population is in the |±4ℏkL ⟩ subspace, weonly keep track of the operation in this subspace. Because of this, we reduce the input state to themodulus square of the matrix elements for only{|2nℏkL ⟩ ⟨±4ℏkL | |n = −3, −2, −1, 0, 1, 2, 3}.

(5.15)The state is chosen so that it forms an experimental observable. This state can be obtained byinitializing the quantum state to one of the |±4ℏkL ⟩ states, applying the control and measuring the
populations for all the momentum states through time-of-flight imaging.
As discussed in the main text, we found that the frequency of the shaking function is important for the performance of the mirror. Simply by modulating at the characteristic frequency,
12ωr , which corresponds to the energy difference between |4ℏkL ⟩ and |2ℏkL ⟩, with a fixed amplitude for a duration of ∼ 3ωr−1 , the channel fidelity can reach as high as 0.8. We take advantageof this knowledge, and try to improve the fidelity on the fixed-amplitude modulation. We choose
our actions to be the amplitude of the sinusoidal modulation to the phase, so the resulting phaseis ϕ(t) = Amp(t) × sin(12ωr t). The amplitudes are chosen from a discrete set of values,
Amp = {0.4, 0.6, 0.8, 1.0, 1.2}.

(5.16)The time interval between each decision point is π/12 ≈ 0.26 ωr−1 , and in each interval the amplitudeis held constant.
The reward function is the same as Eq. (5.9), except that the fidelity F is now defined withthe channel fidelity in the relevant subspace [see Eq. (5.13)].
With the states, actions, rewards specifically designed for reflection, we use the double deep
Q-learning algorithm to learn the strategy for controlling the amplitude of the sinusoidal modulation
to the phase. The deep Q-network used for learning the mirror has an input size of 14, and an
output size of 5, as a result of (5.15) and (5.16). The hyperparameters are listed in Table 5.2.

87
Hyperparameters

Valuesγταϵ decay
episodesHidden size
Batch size
Optimizer

0.99
0.99
0.001
0.0005
8,000
128
32
Adam

Table 5.2: Choice of hyperparameters for the reflection task.

5.4.5

Matter-wave interferometer

We can now cascade together the learned components to make a shaken lattice interferometer.
In between the components, we allow the system to propagate freely with the lattice present butwith ϕ(t) = 0. This gives a time evolution of the wave function in real space as shown in Fig. 5.6,
where the state has been initialized to be in the ground state of the lattice multiplied in real spaceby a Gaussian envelope with a width of a few lattice sites. The free-propagation time is set to
10 ωr−1 . During free propagation, we can see the wave packets propagating at the anticipated±4ℏkL /m group velocity. The anticipated features are evident in the observed splitting, reflection
and recombination of the matter wave.

5.5

Statistical analysis
We compute the final momentum distribution on a grid of acceleration values in order to see

how useful our device is for inertial sensing.
To interpret the outcome of the interferometer, we use Bayes theorem to derive information
about the acceleration from the momentum distribution. This approach is always more effective
than curve-fitting when there is knowledge of the probability distribution generator, and the difference is particularly notable in situations like this when the distributions are multimodal. The
probability distribution of the acceleration after measurement of a particle in a particular momen-

88tum state, P (a|p), is given by,
P (a|p) ∝ P (p|a)P (a) ,(5.17)

where P (a) is the prior distribution of the acceleration, and the proportionality is resolved bynormalizing the total probability to unity. The probability for measuring a particular momentum
p conditioned on the acceleration a, P (p|a) is directly calculated by propagating the wave function
in time. We combine multiple measurements by iterating Eq. (5.17) to formulate the conditional
probability distribution based on the entire measurement record [35]. For a measurement record{p1 , p2 , . . . , pN }, where N is the number of measurements, the probability distribution becomes
P (a|p1 , . . . , pN ) ∝ P (pN |a) . . . P (p1 |a)P (a) .

(5.18)

normalized to unity. Since we assume that each atom is independent of each other, N is thetotal number of atoms we observe. We combine the distributions for each atom using Eq. (5.18)
to accumulate the distribution for a conditioned on the measurement record, and estimate theacceleration by taking the expectation value.
In Fig. 5.7(a), we show an example of the parameter estimation process, where the actual
acceleration is −3 × 10−4 ωr vr (vr = ℏkL /m is the recoil velocity). The corresponding momentum
probability distribution for each atom is shown in the inset. What is important about this distribution is that it is almost unique for each acceleration value and therefore acts as a fingerprint. In
Figure 5.7(b), we show the standard deviation of the acceleration as a function of the number ofmeasurements. A typical number of atoms for an ultracold gas experiment may be of order ∼ 104 ,and with 104 measurements, the estimated acceleration agrees extremely well with the actual value.
We observe that the standard deviation approaches the standard quantum limit, which scales as√∼ 1/ N . This result is consistent with the fact that we are not considering atom interactions inthe system, so each atom may be thought of as making an independent measurement.
We demonstrate that reinforcement learning provides additional capability over traditional
experimental techniques by comparing the sensitivity of the shaken lattice interferometer with theconventional Bragg interferometer. The Bragg interferometry sequence is shown in Fig. 5.1(b).

89
We use the Bayesian parameter estimation result with the Bragg interferometer as a baseline forbenchmarking the shaken lattice interferometer. To make a fair comparison, we consider the casewhere the Bragg interferometry sequence takes the same time as the shaken lattice interferometer
The protocol found by our reinforcement learning algorithm generates a momentum splitting of
8ℏkL (between −4ℏkL and +4ℏkL ) at the beam splitter, four times larger than the 2ℏkL (between−1ℏkL and +1ℏkL ) splitting from Bragg diffraction, and therefore this results in higher sensitivity.
One may argue that it is possible to generate larger momentum splitting with Bragg diffraction,
and while this is true, higher-order transitions typically require much more time than than the shortpulses applied here. For these reasons, we have observed from our simulations that the standard
deviation of the shaken lattice interferometer can be approximately four times lower than that of
the corresponding Bragg interferometer given the time constraint, implying that we have realized
an approximate factor of four in sensitivity gain. The results are shown in Fig. 5.7(b).
We verify that the standard deviations σa for both the shaken lattice interferometer and the
Bragg interferometer are close to the limits set by the Cramér-Rao bound, which can be calculated
from the classical Fisher information as,p
N I1 (a),X 1  ∂P (p|a) 2
I1 (a) =
,
P
(p|a)
∂a
p
σa ≥ 1/(5.19)

where I1 (a) is the Fisher information for one independent measurement. The Cramér-Rao boundfor the shaken lattice interferometer is determined by numerically calculating the classical Fisher
information using the pre-calculated probability P (p|a), and the Cramér-Rao bound for the Bragg
interferometer is determined by the analytical form σa ≥ (2kL T 2 )−1 N −1/2 , where T is the freepropagation time.
Note that the sensitivity presented here is not the ultimate achievable sensitivity since weare constraining system parameters. It would be possible to improve the sensitivity of our interferometer, for example, by increasing the free-propagation time, or by further accelerating thematter-wave components in opposite directions after applying our splitting protocol, as discussed

90
in Refs. [26, 69].

5.6

Conclusion
In summary, we have demonstrated a machine learning methodology to control a complex

quantum system for the purpose of performing a quantum metrology task. Specifically, we demonstrated how to utilize reinforcement learning for the design of a lattice-based matter-wave interferometer. We showed that by shaking the lattice with protocols derived from deep Q-learning,
matter-wave analogs of optical components are realized. We showed that these can be concatenated
together to build a high-precision interferometric device. The multipath interferometer that wasconstructed in this way is capable of the measurement of acceleration with higher sensitivity thanthat achieved by a conventional matter-wave interferometer.
While we have assumed that atom interactions are negligible, they could potentially bringmore quantum advantage if harnessed appropriately. For example, if we could discover a protocol
for shaking the lattice to generate momentum-squeezed states, the sensitivity could potentially behigher than the standard quantum limit that is generated by the statistical averaging of independent atoms. In fact, the quantum entanglement in such squeezed states could allow the resulting
interferometer to approach the ultimate limit to sensitivity set by the Heisenberg uncertainty relation between number and phase. We have focused on demonstrating a specific solution to the
problem of designing a Mach-Zehnder interferometer, but our main outcome is actually to illustrate
an effective design methodology. In the future, this approach could potentially lead to the design ofalternate protocols that construct completely new types of metrological devices with a performancelevel surpassing conventional experimental methods.

91

Figure 5.4: Shaking function for splitting. The learned lattice phase, ϕ(t), (top) is allowed totake on one of a discrete set of five possible values that span the range shown. The momentum
probability distribution (bottom) is initialized to the ground state of the lattice and at the terminal
time well approximates the desired superposition for a beam splitter.

92

Figure 5.5: Shaking function for reflection. The learned lattice phase, ϕ(t), (top) is sinusoidal andthe amplitude is allowed to take on any one of a discrete set of five values at each half-cycle. The
momentum distribution (bottom) is prepared in −4ℏkL and well approximates 4ℏkL at the terminal
time.

93

Figure 5.6: Time evolution of the matter-wave density throughout the entire interferometry sequence. The white dotted lines separate the plot into five regions. In region I, we apply the
splitting protocol. In region II, we allow the matter-wave to propagate freely in the lattice. The
appearance of two wave packets traveling in opposite directions shows that the beam splitter operates as expected. In region III, we apply the mirror shaking function. The matter-wave undergoes
free propagation in region IV again, and the two wave packets switch directions, demonstrating thefunctionality of the mirror. Last, we apply the recombining protocol in region V. Apart from themain closed diamond-shaped paths, we observe that there are auxiliary paths that are fainter but
still clearly evident. They arise due to the imperfection of the components, and also due to theside peaks arising from the third-excited Bloch state.

94

SLI
Bragg
CR boundFigure 5.7: (a) Posterior probability density of the acceleration for the first 100 atoms. The
true acceleration that we aim to reveal by the measurements is −3 × 10−3 ωr vr (red line), andthe measurements are sampled from the momentum distribution at the end of the interferometry
sequence, as shown in the inset. (b) Standard deviation of the acceleration estimated using Bayes
theorem for up to 104 atoms. We show the results from both the shaken lattice interferometer (SLI)
and the Bragg interferometer and conclude that the shaken lattice interferometer has a higher
sensitivity. The standard deviations σa are roughly inversely proportional to the square-root of thenumber of measurements N , for N >√102 . The the black dashed lines are the Cramér-Rao (CR)lower bounds and scale exactly as 1/ N .

Chapter 6

Practical considerations for the shaken lattice accelerometerIn the previous chapter, we demonstrated the construction of a matterwave accelerometerwith individual components designed through reinforcement learning. The sensitivity of the accelerometer shown was obtained from a relatively small enclosed area, with the matterwaves propagating over only about 25 lattice sites. While this was sufficient as a proof of principle demonstration
to show the effectiveness of the matterwave-optics components, there is still a lot of potential togain in the performance of the shaken lattice accelerometer by increasing the dimensions. In fact,
the experiment built by the Anderson group and ColdQuanta Inc. can support an enclosed areathat is of order 100 times larger. On the other hand, we have not taken into account various noise
sources that are unavoidable in a real experiment, which could potentially have detrimental effects
on the performance of the interferometer. In this chapter, we will explore the ultimate limits of the
shaken lattice accelerometer in terms of sensitivity and dynamic range. Furthermore, we will take
into consideration limitations on the engineering design parameters, including laser intensity noise,laser phase noise, dark counts on the CCD camera, and the finite quantum efficiency of the CCD
camera. These simulations provide a guide as to the requirements for equipment specifications toachieve a desired sensitivity.

6.1

From natural units to experimental units
To connect our theory to real-world scenarios, the first step is to convert the natural units usedin the theory to experimental units. The system we consider is made up of an optical lattice formed

96
by counter-propagating lasers with wavelength λL = 852 nm, and the atoms trapped in the lattice
are rubidium atoms 87 Rb, whose individual masses are m = 87 amu. With ℏ = 1.05×10−34 m2 ·kg/s
and 1 amu = 1.66 × 10−27 kg, we can calculate the recoil frequency ωr and the recoil velocity vr ,
ωr =
vr =

ℏkL2
= 1.98 × 104 s−1 ,2mℏkL= 0.0054 m/s ,
m

(6.1)
(6.2)

where we have used kL = 2π/λL . Note that the recoil frequency is expressed in radians per second
above and converts to 3.15 kHz in hertz. The natural unit we have used for acceleration is ωr vr ,and from Eqs. (6.1-6.2), we get
1 ωr vr = 106 m/s2 ≈ 10 g ,(6.3)

where g is the gravitational acceleration on the surface of Earth. Time has been presented in ωr−1 ,where 1 ωr−1 = 50 µs, and the unit for spatial dimension, kL−1 , converts to 1 kL−1 = 136 nm.
In Chapter. 5, the sensitivity for the acceleration with 104 measurements in the shaken latticeinterferometer reached 10−5 ωr vr = 1.1 × 10−3 m/s2 ≈ 10−4 g. The total time for one run of theinterferometry sequence was Ttot ≈ 27 ωr−1 = 1.36 ms, during which the wave function traveled
across a distance of about L = 160 kL−1 = 0.0217 mm.

6.2

Noise sources in the experiment

6.2.1

Laser linewidth

A laser with non-zero linewidth has phase diffusion that can be modeled as white noise. Thismeans the phase undergoes a random walk in time, and the phase noise can be described by a
Wiener process,
ϕt+dt − ϕt ∼ N (0, ∆ dt) ,(6.4)

where N (µ, σ 2 ) describes a Gaussian distribution with mean µ and standard deviation σ, ∆ is thelaser linewidth, and dt is the time step used in the simulation. If the two counter-propagating lasers
have common phase noises, which is the case when the counter-propagating light is generated from97
a common source, e.g., by a single retro-reflected laser, then the relative phase between them, i.e.,
ϕ in Eq. (5.1), will not suffer from the noise from the laser linewidth. However, if the laser noises
are independent, which is the case when the lattice is formed by two separate light sources, thenwe need to take the phase noise into account. The standard deviation of the noise in the relative√
phase will be 2 times the standard deviation in a single laser.
To incorporate the laser noise in our simulations, we assume thatϕ(t) = ϕC (t) + ϕr (t) ,(6.5)

where ϕC is the controlled phase, and ϕr is the phase noise that fluctuates as described in Eq. (6.4).
The controlled phase functions are the ones the reinforcement learning agent learns in a noiselessenvironment, as shown in Chapter 5. We simulate the output momentum distributions from the interferometer with ∆ = 2×10−4 ωr , which corresponds to a linewidth of ∼ 16Hz, and use a Bayesian
approach to reconstruct the acceleration. The value of ∆ is chosen to give the same accumulated
phase during the interferometer sequence that would be generated by the smallest detectable acceleration in the simulation with monochromatic lasers. Fig. (6.1)(a) shows the reconstruction process
for the acceleration in this case. It may be observed that even though the distribution becomesnarrower as the number of measurements increases, the center of the peak is shifted from the true
value, meaning that even though the precision converges, the estimation is biased, compromising
accuracy. This can also be concluded from Fig. (6.1)(b) and (c).
Fig. (6.1)(b) shows the standard deviation of the acceleration, which goes down as the standard quantum limit. Fig. (6.1)(c) shows the mean squared error of the estimated acceleration. The
value where σa saturates corresponds to the bias of the estimation.
The non-zero laser linewidth also results in a spread in the wave vector kL . However, since∆ν/ν ∼ ∆kL /kL ∼ 3 MHz/300 THz ∼ 10−8 , it is expected that the fluctuation in the wave vector
will not be a main source of error.

98

Figure 6.1: (a) Bayesian reconstruction of the true acceleration in the presence of a phase noisethat corresponds to a linewidth ∆ = 2 × 10−4 ωr . The red line indicates the true acceleration.
(b) Standard deviation of the reconstructed acceleration. (c) Mean squared error of the acceleration.

99
6.2.2

Laser intensity fluctuation

To model the laser intensity noise, we assume that the lattice depth fluctuates as
V0 (t) = V0 [1 + W (t)] ,(6.6)

where V0 represents the average lattice depth, and W (t) is a white noise that fluctuates as
Wt+dt − Wt ∼ N (0, ξ dt) .

(6.7)

We simulate the shaken lattice interferometer with laser intensity fluctuation rate ξ = 2 × 10−4 ωr ,and reconstruct the acceleration signal using a Bayesian approach. This rate can be interpreted as
an empirical upper bound, and can be expected to be much larger than what would be typical inexperiment. Fig. (6.2)(a) shows the reconstruction process for the first 100 measurements (atoms).
Fig. (6.2)(b) shows the standard deviation, which still scales as the standard quantum limit atlarge N despite the large level of noise at small number of measurements. However, in Fig. (6.2)(c),we can see that the level of the mean squared error saturates at a relatively large value. The
explanation is that during the first 100 measurements, there the distribution of the acceleration
is multi-peaked due to shot noise. After we make more measurements, the distribution becomessingle-peaked, which is a result of the maximum likelihood principle, and therefore, the standard
deviation starts to decrease. However, the peak may not coincide with the true value due to thepresence of the intensity noise, and therefore, we observe the high level of mean squared error.
In addition to the white noise that is uniform for all frequencies, we also simulate the intensity
noise with a certain pattern in the power spectral density, which can be captured on a photodetector
connected to a spectrum analyzer. The resulting time series is shown in Fig. 6.3(a). We only use
a segment of the time series for each acceleration value, and for the next value we use the next
segment, and so on. When we reach the end of the series, we go back to the beginning. This ensuresa certain amount of randomness. The momentum distribution of the interferometer is shown in
Fig. 6.3(a). The maximum difference between the noisy distribution and the noiseless distribution
is only about 1%.

100

Figure 6.2: (a) Reconstruction of the true acceleration in the presence of intensity noise ξ =2 × 10−4 ωr . (b) Standard deviation of the acceleration. (c) Mean squared error of the acceleration.

101

(c)

Figure 6.3: (a) Reconstruction of the acceleration from a system with laser intensity noise as shownin (c). (b) The standard deviation of the acceleration as calculated from the reconstruction. (c)Time series of the noise. The noise is defined as the ratio between the intensity fluctuation and themean intensity.

102
6.2.3

Quantum efficiency

An important property of the photodetector is that it may not detect all the atoms. Thiseffect is known as the quantum efficiency and is considered to be an imaging noise. Suppose
the detector has a quantum efficiency of E, then each atom has a probability of 1 − E to be notcounted. We simulate this effect by ignoring each sample drawn from the momentum distribution
with probability 1 − E. The result of this is summarized in Fig. (6.4).
6.2.4

Dark counts

Dark counts come from the background noise on the CCD image when there is actually no
atom. We simulate this effect by adding uniform random noise on top of the original samples. The
noise level is defined by the ratio between the number of random samples and the number of atoms
there actually is. The result of this is summarized in Fig. (6.5).

6.3

Dynamic range
To study the limit for the dynamic range, we simulate the interferometry outcome for a wide

range of acceleration. We observe that within a range of acceleration close to 1 g, the momentum
distribution for each acceleration value is still unique, and therefore acts as a fingerprint that allows
us to make unbiased estimations. The results are shown in Fig. 6.6.
We then estimate the acceleration from the interferometry outcome simulated using the true
acceleration. The estimation is calculated by taking the mean from the reconstructed probability
distribution of the acceleration. We do this for all the acceleration values in the range presented
in Fig. 6.6, and look at the relation between the true acceleration and the estimated acceleration.
The fact that the curve lies on x = y means that the estimation is unbiased in this range.

103

Figure 6.4: (a) Reconstruction of acceleration with E = 0.9. (b) Standard deviation. We assume
that there are originally 10000 atoms, but since E ≠ 1, we only get 9035 counts on the detector.
The inset shows a zoom-in of N = 8000 to N = 9035.

Figure 6.5: Error and standard deviation in the presence of dark counts. The absolute value of theerror grows larger as the noise level increases, whereas the standard deviation remains almost the
same.

104

Figure 6.6: Momentum distribution for a wide range of acceleration. The free-propagation time isset to 10 ωr−1 .

105
(a)

(b)

(c)

Figure 6.7: (a) A plot of the estimated acceleration versus the true acceleration. Each accelerationis reconstructed from 10000 samples. The linear relation verifies that the estimations are unbiased.
(b) A scatter plot of the error. It can be seen that the error is uniformly distributed across thewhole range of true acceleration. (c) A histogram of the error, where the width of each bin is
4 × 10−6 (ωr vr ).

106

6.4

The ultimate limit for sensitivity
Now let us extend the free propagation time to T = 1000 ωr−1 . This corresponds to a freepropagation distance of ∼ 1 mm, which fits into a typical optical lattice experiment in a glass cell.
We found that at this limit, the sensitivity reaches 10−9 ωr vr = 10−8 g, and if we repeat theinterferometry 100 times, which takes about 20 min due to the finite loading time for the BEC, we
could push the sensitivity to 10−9 g. At this level, there are already some interesting applications forspace navigation, such as detecting the sudden change in acceleration caused by the solar radiationpressure on a spacecraft when it emerges from a region of shadow.

6.5

Outlook
Since the noises could impose a bias to the parameter estimation as shown in our analysis, it is

crucial to resolve this problem, and machine learning could potentially mitigate this adverse effect.
One direction is to employ a supervised learning algorithm that learns the relation between the
momentum distribution and the acceleration from the noisy data. The data are labelled by the true
acceleration, which is considered a known parameter under a well-controlled experimental setup,
and after the supervised-learning model is trained, we could test its efficacy with distributions thatcome from an unknown acceleration. The other direction would be to develop a reinforcement
learning agent that is capable of learning in a noisy environment and finding the protocols thatare the most robust against the noises. In this case, we could add noise robustness to the reward
function, so that the reward is large for small variations of the distribution with respect to laser
intensity or frequency. The learning cycle would then result in the agent optimizing with respect
to this quantity.

107
(a)

(b)

(c)

(d)

Figure 6.8: (a) Final momentum distribution for an interferometer where the free propagation timeis 1000 ωr−1 . (b) A zoom-in of (a) near a = 0. (c) The error and the standard deviation of theestimation of the acceleration from plot (b) with 10000 samples. (d) Histogram of the estimation
error.

Chapter 7

How to train your gyro: Reinforcement learning for rotation sensing with ashaken optical Lattice

7.1

Introduction
If one asks an expert for advice about how to measure tiny displacements or forces over large

distances, they will likely point to the very large body of literature on interferometry, perhaps
highlighting LIGO [1] and other such remarkable interferometry systems. Optical interferometers
are not the only ones useful in high precision measurements, and indeed matter-wave interferometers
have also offered the potential for exquisite sensitivity to inertial forces. In the work by Weidnerand Anderson, a new kind of matter-wave interferometric accelerometer that utilizes an opticallattice was proposed and demonstrated experimentally [102, 100, 101]. The purpose of the opticallattice is to provide robustness of the system in the face of a harsh dynamical environment typicalof real-world applications. It has been demonstrated that quantum design of optical-lattice-based
interferometers can benefit greatly from machine learning approaches, since the machine can explore
an extensive number of possible control protocols and come up with one that is not intuitive butvery effective. The previous chapter established that one can ‘teach’ the system to become sensitiveto acceleration by learning how to modulate the phase of the lattice with unconventional patterns,
that is, how to shake the lattice such that the system carries out an interferometric protocol.
While this work demonstrated the power of machine learning to construct an interferometer tosurpass conventional Bragg interferometry, the machine learning was constrained, let us say, by
an expert’s understanding of how interferometry is to take place. Indeed, nearly every kind of

109
interferometry, optical or matter wave, takes place as a sequence of wavefront splitting, reflection,and recombination, with wave propagation in-between, and this is known as the Mach-Zehnderconfiguration [54, 107].
If we free the learning agent from the constraints set by conventional configurations, it could
potentially learn beyond what human has explored. For example, Ref. [57] has shown that an
active-learning agent can predict novel configurations for quantum optics experiments that generatehighly-entangled states. This chapter presents results from the application of machine learning totrain a two-dimensional optical lattice to sense rotation, i.e., to form a gyroscope, as depicted in
Figure 7.1. In contrast to the accelerometer, in this case we utilize reinforcement learning in amanner that removes the prejudice of an expert about how interferometry is supposed to be done.
Instead, learning is driven by a focus on performance metrics — that is, what it means to be a good
gyroscope rather than what it means to be a proper interferometer. We show that the solution
looks nothing like standard interferometry and that the system, in fact, exhibits a sensitivity torotation that is multiple factors better than a standard interferometer. The fact that there is anumerical improvement per se is perhaps less important than the fact that the lattice interferometer
exhibits the complexity enabled by quantum mechanics; one associated with a large Hilbert spacedimensionality. The learned solution outperforms the expert’s preconceived notions about how one
might use that complexity to advantage. Moreover, the machine learning can cater to a set of
performance metrics, such as a combination of sensitivity, dynamic range, and tolerance to noiseand experimental drifts. The learning algorithm navigates a trade-space that may well be opaque
to the expert designer.

7.2

Physical model
The system we consider is one where the atoms are loaded into a two-dimensional opticallattice. The lattices in both directions can be ‘shaken’ by varying the phases of the corresponding
interfering laser beams. The entire system, shown in Figure 7.1, is considered in a rotating non-

110

e#
Shak

e"
ak
Sh

!
Figure 7.1: A two-dimensional optical lattice potential in a rotating frame, where atoms are trapped
and their dynamics can be controlled by the shaking of the lattice. The blue arrows indicate possible
trajectories associated with the splitting of the matterwaves.

111
inertial frame, and therefore described by the Hamiltonian
Ĥ =

Vy
p̂2
Vx
− Ω · L̂ −cos(2kL x̂ + ϕx ) −cos(2kL ŷ + ϕy ) ,
2m
2
2

(7.1)

where m is the mass of a single atom, kL is the laser wavenumber, ϕx and ϕy are the phase differences
between the two counter-propagating lasers in the x and y directions respectively, and Vx and Vy
are the corresponding lattice potential strengths. The term Ω · L̂ describes the rotational kinetic
energy of the system, with L̂ = x̂ × p̂ the angular momentum, and Ω being the angular velocity.
It is the magnitude, Ω, that we aim to measure with high precision.

7.3

Two-path Sagnac matter-wave interferometer
To measure the angular velocity, we first consider a two-path loop configuration [32, 50]

that resembles the construction of the fiber-optic gyroscope [93], which directly utilizes the Sagnac
effect [77, 6]. The Sagnac effect describes the physical process that generates a phase difference
between the two waves that propagate in the opposite directions in a loop when in a non-inertial
rotating frame of reference. In general, this predicts the property that the phase difference, ∆ϕ,is proportional to both the angular velocity and to the area that the two paths enclose, A. In thecase of matter-waves, this relationship can be derived to be ∆ϕ = 4mΩ · A/ℏ.
A two-path Sagnac matter-wave interferometer can be constructed in the lattice using thebeam-splitting and reflecting protocols presented in [17], based on components optimized by reinforcement learning. In that case, the splitting protocol was demonstrated to transfer the groundstate to a superposition of the |±4ℏkL ⟩ states, and the reflecting protocol to map any linear combination of the |±4ℏkL ⟩ states to the corresponding combination of the |∓4ℏkL ⟩ states.
In order to build a gyroscope out of these components, the sequence of operations is thefollowing. Initially, the wavefunction is prepared in the ground state of the 2D lattice. In thex-direction, the 1D beam-splitting protocol is applied, allowing free propagation for a duration,
denoted as T , and then the reflecting protocol is applied. Following this, free propagation occurs
for a duration 2T , and then the reflecting protocol is applied, followed by free propagation for112a duration T , and finally the recombining protocol. In the y-direction, the lattice operates as aconveyor belt for the trapped atoms. The lattice is smoothly accelerated to the velocity 4ℏkL /maccording to the adiabaticity criterion, and then translated at constant velocity. During the halfway
point of the sequence, when the two-paths in the x-directions cross, the velocity of the y-lattice issymmetrically decelerated through zero to −4ℏkL /m, so that the lattice can be translated backward
at the constant velocity from that point on. The final step is to adiabatically accelerate the y-lattice
back to zero velocity for the final recombination protocol.
To ensure that the acceleration of the y-lattice is indeed adiabatic, let us investigate theadiabaticity condition for the wavefunction to stay in the ground state of an accelerating lattice.
The Hamiltonian in the y direction in the lab frame is
Ĥy =


p̂y
Vy
−
cos 2kL ŷ + ϕ(t) .
2m
2

(7.2)


With the unitary transform by the operator Û = exp − ip̂y ϕ/2ℏkL , we move into the frame
co-moving with the lattice. The Hamiltonian in the co-moving frame is

p̂y
Vy
ϕ̇
ˆe
H
−
cos 2kL ŷ + ϕ(t) +
p̂y .
(7.3)
y =
2m
2
2kL
E



The instantaneous eigenstates are ψeN,q = exp imŷ ϕ̇/2ℏkL exp − ip̂ϕ/2ℏkL |ψN,q ⟩, where|ψN,q ⟩ are the Bloch eigenstates with band index N and quasimomentum q. If we wish the wavefunction stays in the ground state (N = 0), then the quantum adiabatic theorem requires that [19]X ℏ ⟨m| p̂ ϕ̈/2kL |0⟩
X ℏ ⟨m| He˙ |0⟩=≪ 1.
2(Em − En )(Em − En )2

m̸=0

(7.4)

m̸=0

For Vy = 10 Er , numerical calculations show that the adiabaticity condition is |ϕ̈| ≪ 60 ωr2 . In our
simulations, we use |ϕ̈| = 6.4ωr2 , which satisfies the adiabaticity condition.
If higher sensitivity is desired, one could trivially extend this sequence to allow the matterwaves to travel around in multiple loops before applying the last recombination step, although inpractice, one needs to balance this against the resulting reduced fringe visibility. A schematic to
illustrate the paths that the wave functions traverse is shown in Fig. 7.2. The point here is to113

Figure 7.2: The clockwise path (blue) and the counter-clockwise path (red) taken by the wave
functions in the two-path shaken lattice gyroscope. The vertical back line is where the wave
functions were initially split and are eventually recombined.

114
demonstrate that a gyroscope can indeed be constructed in a two-dimensional lattice, and to set
a baseline for sensitivity for the purpose of making comparison with the gyroscope designed by
end-to-end reinforcement learning that will follow later.
The sensitivity of a gyroscope constructed this way is presented in terms of the classicalFisher information in Fig. 7.3. The classical Fisher information is a quantitative measure of thesensitivity of the output distribution on the metrological parameter. Note that the value of theFisher information presented here is normalized to unity for an ideal two-path gyroscope with thesame momentum splitting and free-propagation duration as the one we construct, i.e., (4p0 T 2 )2 ,where p0 = 4ℏkL . The shaken lattice gyroscope gives a Fisher information of order unity as can be
anticipated, and the slight deviation is due to imperfections in the components and the fact thatthe shaken lattice is does not fit perfectly within the two-path framework.

7.4

Effective 1D model
One challenge of simulating the two-dimensional rotating lattice is that, in general, the com-putational complexity scales quadratically with the number of momentum states. To enable efficientsimulations, which is especially important when running many machine-learning epochs, we reduce
the system to a linear complexity by introducing an approximation of a deep potential in they-direction that reduces the system to an effective one-dimensional model. Explicitly, we take advantage of the fact that when the deep lattice in the y-direction is adiabatically translated, theshape of the wavepacket in the y-direction transported along with the lattice remains approximately invariant. This approximation is encapsulated by a separable ansatz, where we assume thetwo-dimensional wavefunction is a simple product state as represented by the tensor product of the
wavefunctions in the x and y directions, i.e.,

115

Figure 7.3: The Fisher information for the machine-learned shaken-lattice gyroscope normalizedwith respect to the Fisher information that results from the gyroscope in which the matterwaves
follow the ideal paths of the conventional two-path inteferometer.

116|ψ⟩ ≈ |ψx ⟩ ⊗ |ψy ⟩ ,p̂y ≈ ⟨py ⟩ ,ŷ ≈ ⟨y⟩ .

(7.5)

With this ansatz, the time-dependent solution and the interference that is formed is determined by
the x-direction evolution, modeled from,
ex =
H

V0
p̂x
ez ,
−
cos(2kL x̂ + ϕx ) − Ω L
2m
2

(7.6)


e z = ⟨py ⟩x̂ − ⟨y⟩p̂x . The mean position and momentum in the y-direction for the sequence
where Ldescribed in Sec. 7.3 is shown in Fig. (7.4). The oscillations in the momentum plot are due to theatoms being trapped in the bottom of the lattice sites, which are close to harmonic traps, and theaccleration phase being not perfectly adiabatic.
Let us further discuss the validity of the separation ansatz. We can write the 2D wavefunctionin the form of its Schmidt decomposition,
ψ(x, y) =Xqp(i) ψx(i) (x) ψy(i) (y)

(7.7)

i

This decomposition can be obtained numerically from singular value decomposition. If the highest singular value p(0) is close to identity, we can drop the other terms in the summation, andapproximate the 2D wavefunction as the separable form, i.e.,ψ(x, y) ≈ ψx(0) (x) ψy(0) (y)(7.8)

To verify that in the case where the y-wavepacket is adiabatically transported along the lattice, theapproximation shown above is valid, we compute the Schmidt decomposition of the 2D wavefunction evolved from the full 2D Hamiltonian. The highest singular value is 0.95. We compare the´population distribution in x from the full 2D wavefunction, |ψ2D (x)|2 = |ψ(x, y)|2 dy, and that(0)from the Schmidt decomposition, |ψx (x)|2 , and as shown in Fig. 7.5, the two populations agree
well.

117

Figure 7.4: (a) Average position of the wavepacket in the y-direction. (b) Average momentum of
the wavepacket in the y-direction. The lattice depth Vy is set to 10 ωr .

Figure 7.5: Comparison between the simulations results of the probability distributions in thex direction from the effective 1D model (blue), the full 2D model (orange), and the highest singularvalue component from the 2D model (green).

118
The validity of the approximation allows us to use the separable ansatz, i.e.,ψsep (x, y) = ψx (x)ψy (y)(7.9)

and evolve the wavefunctions in the x and y directions separately using the effective 1D Hamiltonians. From now on, we will work under the effective 1D assumption, where only the dynamics in
x-direction is investigated.

7.5

End-to-end design with reinforcement learning
We have shown how to build a two-path Sagnac gyroscope from individual interferometercomponents. In this section, we consider applying reinforcement learning to the design of thegyroscope to go beyond the traditional Mach-Zehnder configuration. The details of the principles
of reinforcement learning and deep-Q learning can be found in [17] and will not be repeated here.
Our purpose is to demonstrate how to formulate an end-to-end design problem for a matter-wave
gyroscope within the scope of the reinforcement learning framework.
A basic reinforcement learning framework consists of a agent that select actions based on theobserved state of the environment, and an environment that interacts with the action and generates
a reward based on the observed quality of the current state. In our design framework, the state(not the quantum state, but the state of the environment within the context of machine learning)
contains information about the population in the discrete momentum basis in x and the average
position and momentum in y. The actions are defined as a finite set of discrete choices for the time
derivative of the phase, i.e., ϕ̇. This choice is motivated by the experimental implementation forthe control of the lattice phase, where an acoustic-optical modulator is employed and the frequencyof the modulation is proportional to ϕ̇. Let us now elaborate more on the states and actions wehave chosen. The state given to the reinforcement learning agent consists of information about the
populations in {|2nℏkL ⟩ | ± 3, ±2, ±1, 0} in the x-direction and the position and momentum in they-direction. It is represented by the vector,⟨ŷ⟩ ⟨p̂y ⟩s = P−3 , P−2 , P−1 , P−0 , P−1 , P−2 , P−3 ,
,
,
40 4

(7.10)

119
where Pn represent the population in the |2nℏkL ⟩ state. The action is the time derivative of thelaser phase in the x-direction selected from a discrete set of values,A ∈ {−ϕ̇ = −8, ϕ̇ = 0, ϕ̇ = 8}.

(7.11)In order to compute the reward, we construct the environment as three copies of the quantum
system with different rotation rates, Ω = [Ω0 − dΩ, Ω0 , Ω0 + dΩ], for small dΩ and with each
evolving according to Eq. (7.6). This choice is motivated by the requirement to compute thederivative with respect to Ω by numerical symmetric finite differencing. Explicitly, the reward isdefined as the classical Fisher information computed with respect to Ω using the following formula,I1 (Ω) =

X
p



1
∂Pr(p|Ω) 2.
Pr(p|Ω)∂Ω(7.12)The goal of learning is to find a sequence of ϕ̇ that optimizes the Fisher information at the terminaltime.
We point out that the Fisher information is related to the theoretical bound (Cramér-Rao
bound) for the standard deviation of the measurement for Ω as follows,pσΩ ≥ 1/ N I1 (Ω),(7.13)

so by optimizing the Fisher information with reinforcement learning, we are effectively optimizing
the measurement sensitivity with respect to the generated quantum state at the unitary output,
but based only on a specific choice of the angular frequency, Ω0 . A variety of extensions to what
we present here are clearly possible, based on alternative design conditions.

7.6

Results for end-to-end design
Within the context of the design framework described in Sec. 7.5, we are able to obtain lattice

shaking protocols that outperform the conventional two-path gyroscope demonstrated in Sec. 7.3 interms of sensitivity. The protocol varies for every instantiation due to the randomized initializationof the agent, and Fig. 7.6 shows one of the machine-learned protocols. The time-evolution of the120
resulting machine-learned multi-path interferometer is shown in Fig. 7.6. The pattern is reminiscentof the speckle patterns that emerge from a multi-mode fiber (fiber specklegram), which are known tobe sensitive interference detectors of inertial phase. On the other hand, in the case of a multi-mode
fiber, that patterns are scrambled by temperature changes or strain on the fiber. For the multimode matter-wave interferometer we consider here, the situation is quite different, and in spite of
the irregular pattern it generates, we show that this device can, in principle, measure rotationalsignals in an unbiased fashion and reach the sensitivity predicted by the theoretical bound.
The momentum distributions that are produced by the multi-path interferometer under conditions of different values of the rotation rate are shown in Fig. 7.7. The larger the variation in the
momentum distribution is as the rotation rate is varied, the more sensitive is the interferometer.
The Fisher information calculated from the distribution is shown in Fig. 7.7. The value on the
y-axis is the ratio between the Fisher information from the multi-mode interferometer constructedby reinforcement learning and the one shown in in Sec. 7.3. We optimize the Fisher information
around Ω0 = 0, and one can see that the Fisher information is maximum at Ω = 0, with the
dynamic range being sacrificed for the high sensitivity. The machine-learned gyroscope achieved aFisher information of around 25 times higher than the two-path gyroscope, meaning that the sensitivity is improved by an extremely impressive factor of 5. If we make a comparison to a gyroscope
constructed from Bragg diffraction, the improvement is 4 × 5 = 20 − fold.
To extract the rotational signal from the distribution shown in Fig. 7.7, we apply Bayesianreconstruction. This means that we iteratively update the prior distribution for Ω from eachsample. The reconstruction of the rotational signal is shown in Fig. 7.8(a), and the estimation forΩ is unbiased. In Fig. 7.8(b), we plot the standard deviation of Ω with the number of samples
taken, and see that the sensitivity scales inversely with the square root of the number of atoms, as
expected from independent measurements.

121

Figure 7.6: Phase and the spatial wavefunction as a function of time. The shaking protocol is
generated from the reinforcement learning agent with ϕ̇ as the actions. There are 125 steps in thisprotocol. The wavefunction is initialized in the ground state with a momentum width σp = 0.1,
and evolve with the shaking protocol shown above.

122

Figure 7.7: Population distribution in the momentum basis and the resulting classical Fisher information for the machine-learned shaken-lattice gyroscope as a function of the rotation rate Ω.
The value of the Fisher information is normalized with the Fisher information of a conventional
two-path interferometer with splitting ∆p = 8ℏkL .

123

RL
CR boundFigure 7.8: (a) A reconstruction of the rotation rate from 1 sample to 100 samples. (b) The
standard deviation of the rotation rate as a function of the atom number.
√ The result at large atom
number is consistent with the Cramér-Rao bound, which scales as 1/ N .

124

7.7

Conclusion and outlook
In this chapter, we have proposed and analyzed the efficacy of the use of deep reinforcement

learning for the optimization of a quantum sensor. We have demonstrated the explicit application ofthis concept through an example of the design of a shaken-lattice matter-wave gyroscope. We first
showed that a Mach-Zehnder-type matter-wave gyroscope can be constructed in a 2D lattice froma conveyor belt in the y-direction and the previously developed beam-splitting and the reflecting
shaking functions in the x-direction. We then applied reinforcement learning in an end-to-end
fashion to find a shaking protocol that improves the sensitivity of the gyroscope by a factor of 5compared to the conventional Mach-Zehnder interferometry. These impressive results open up the
possibility for applying the same principal machinery to a variety of quantum sensing problems
with more complex landscapes, such as multi-parameter estimation and entanglement-enhanced
metrology. However, the implications of our results go beyond this framework, and the same ideasof reinforcement learning may be applied to a variety of quantum and classical systems where designof complex circuits is needed.

Chapter 8

Summary and outlook

In this thesis, we applied the basic concepts of quantum metrology and machine learning,presented in Chapter 2 and 3 respectively, to the design of quantum systems for extreme sensing.
We gave three examples of quantum design in this thesis corresponding to three levels of design philosophy. First, we showed that the periodic modulation of the scattering length in a Bose-Einsteincondensate at the resonant frequency can generate squeezed states in momentum, and proposed a
protocol to generate Ramsey fringes to probe the coherence in this system. Next, we demonstrated
the application of reinforcement learning to the design shaking functions that generate the essential
components of Mach-Zehnder atom interferometry in a one-dimensional lattice. We were able toconstruct an interferometer with a larger enclosed space-time area than the conventional Bragg
interferometer, and therefore achieved a higher sensitivity to acceleration. We further investigated
how well this interferometer would work in a hypothetical experimental setting. Finally, we moved
forward and used reinforcement learning in an end-to-end manner, with the aim to design a gyroscope in a two-dimensional shaken lattice that is not constrained by our prior knowledge about
how interferometers should work. We discovered the shaking function that led to a device that
performs better than the conventional Sagnac gyroscope in terms of sensitivity to rotation.
The reinforcement learning framework for quantum design is general, and has inspired proposals for many potential future research works, including:• Using the formalism established in Chapter 4 to simulate an environment that can be controlled by a reinforcement learning agent, with the goal of exploring methods to generate

126
squeezed states in a Bose-Einstein condensate that can potentially be used for entanglementenhanced atom interferometry.
• Extending the design of a unitary operator demonstrated in Chapter 5 to the design ofa general trace-preserving quantum channel in the context of an open quantum system,
which could be useful for tackling noise in a quantum processor.
• Designing the optimal encoding and decoding protocols for a multi-qubit programmablequantum sensor using natively available gates on platforms such as ion or neutral atomarrays, for the purpose of achieving entanglement-enhanced metrology.
• Finding the protocol to prepare the optimal quantum state for a multi-parameter estimation
scheme, in particular for a set of non-commuting observables.
• Improving the momentum capture range in laser cooling and particle slowing, which leads
to a larger number of particles that can be used for performing a metrology task, andtherefore increasing the sensitivity.
Even though the central focus of this thesis has been on metrology, the design philosophy is
general, and can now be applied to a wide variety of problems relevant to engineering quantum systems to perform desired tasks. These include algorithmic operations that are relevant for quantum
computation, quantum communication applications, quantum simulation ideas and concepts, andquantum control, and many of these contain quantum entanglement and many-body correlations
that will go beyond what we have presented here.

Bibliography[1] B. P. Abbott et al. Observation of gravitational waves from a binary black hole merger. Phys.
Rev. Lett., 116:061102, Feb 2016.
[2] Z. An and D. L. Zhou. Deep reinforcement learning for quantum gate control.
(Europhysics Letters), 126(6):60002, jul 2019.

EPL[3] G. B. Arfken, H. J. Weber, and F. E. Harris. Chapter 14 - bessel functions. In G. B.
Arfken, H. J. Weber, and F. E. Harris, editors, Mathematical Methods for Physicists (Seventh
Edition), pages 643 – 713. Academic Press, Boston, seventh edition edition, 2013.
[4] M. Arratia. On the jets emitted by driven bose–einstein condensates. Journal of Physics B:
Atomic, Molecular and Optical Physics, 52(5):055301, feb 2019.
[5] B. Barrett, A. Bertoldi, and P. Bouyer. Inertial quantum sensors using light and matter.
Physica Scripta, 91(5):053006, apr 2016.
[6] B. Barrett, R. Geiger, I. Dutta, M. Meunier, B. Canuel, A. Gauguet, P. Bouyer, and A. Landragin. The sagnac effect: 20 years of development in matter-wave interferometry. Comptes
Rendus Physique, 15(10):875–883, 2014. The Sagnac effect: 100 years later / L’effet Sagnac: 100 ans après.
[7] R. Bellman. On the theory of dynamic programming. Proceedings of the National Academyof Sciences, 38(8):716–719, 1952.
[8] T. Bergeman, M. G. Moore, and M. Olshanii. Atom-atom scattering under cylindrical harmonic confinement: Numerical and analytic studies of the confinement induced resonance.
Phys. Rev. Lett., 91:163201, Oct 2003.
[9] J.-P. Blaizot and G. Pipka. Quantum theory of finite systems. MIT Press, 1986.
[10] N. N. Bogolyubov. On the theory of superfluidity. J. Phys. (USSR), 11:23–32, 1947.
[11] C. Bordé. Atomic interferometry with internal state labelling. Physics Letters A, 140(1):10– 12, 1989.
[12] T. Bothwell, C. J. Kennedy, A. Aeppli, D. Kedar, J. M. Robinson, E. Oelker, A. Staron, andJ. Ye. Resolving the gravitational redshift across a millimetre-scale atomic sample. Nature,602(7897):420–424, Feb 2022.

128[13] M. Bukov. Reinforcement learning for autonomous preparation of floquet-engineered states:
Inverting the quantum kapitza oscillator. Phys. Rev. B, 98:224305, Dec 2018.
[14] M. Bukov, A. G. R. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta. Reinforcement
learning in different phases of quantum control. Phys. Rev. X, 8:031086, Sep 2018.
[15] J. Carrasquilla. Machine learning for quantum matter. Advances in Physics: X, 5(1):1797528,
2020.
[16] T. Chen and B. Yan. Characterization of stimulated excitations in a driven bose-einsteincondensate. Phys. Rev. A, 98:063615, Dec 2018.
[17] L.-Y. Chih and M. Holland. Reinforcement-learning-based matter-wave interferometer in a
shaken optical lattice. Phys. Rev. Research, 3:033279, Sep 2021.
[18] L. W. Clark, A. Gaj, L. Feng, and C. Chin. Collective emission of matter-wave jets fromdriven bose-einstein condensates. Nature, 551:356 EP –, Nov 2017.
[19] D. Comparat. General conditions for quantum adiabatic evolution. Phys. Rev. A, 80:012106,
Jul 2009.
[20] G. M. D’Ariano and M. G. A. Paris. Arbitrary precision in multipath interferometry. Phys.
Rev. A, 55:2267–2271, Mar 1997.
[21] E. A. Donley, N. R. Claussen, S. L. Cornish, J. L. Roberts, E. A. Cornell, and C. E. Wieman.
Dynamics of collapsing and exploding bose-einstein condensates. Nature, 412(6844):295–299,
2001.
[22] L. Feng, J. Hu, L. W. Clark, and C. Chin. Correlations in high-harmonic generation ofmatter-wave jets revealed by pattern recognition. Science, 363(6426):521–524, 2019.
[23] T. Fösel, P. Tighineanu, T. Weiss, and F. Marquardt. Reinforcement learning with neuralnetworks for quantum feedback. Phys. Rev. X, 8:031084, Sep 2018.
[24] S. Fray, C. A. Diez, T. W. Hänsch, and M. Weitz. Atomic interferometer with amplitude
gratings of light and its applications to atom based tests of the equivalence principle. Phys.
Rev. Lett., 93:240404, Dec 2004.
[25] H. Fu, L. Feng, B. M. Anderson, L. W. Clark, J. Hu, J. W. Andrade, C. Chin, and K. Levin.
Density waves and jet emission asymmetry in bose fireworks. Phys. Rev. Lett., 121:243001,
Dec 2018.
[26] M. Gebbe, J.-N. Siemß, M. Gersemann, H. Müntinga, S. Herrmann, C. Lämmerzahl,
H. Ahlers, N. Gaaloul, C. Schubert, K. Hammerer, S. Abend, and E. M. Rasel. Twin-lattice
atom interferometry. Nature Communications, 12(1):2544, 2021.
[27] D. M. Giltner, R. W. McGowan, and S. A. Lee. Theoretical and experimental study of thebragg scattering of atoms from a standing light wave. Phys. Rev. A, 52:3966–3972, Nov 1995.
[28] S. Giorgini. Collisionless dynamics of dilute bose gases: Role of quantum and thermal fluctuations. Phys. Rev. A, 61:063615, May 2000.

129[29] V. Giovannetti, S. Lloyd, and L. Maccone.
Photonics, 5(4):222–229, Apr 2011.

Advances in quantum metrology.

Nature[30] E. Greplova, C. K. Andersen, and K. Mølmer. Quantum parameter estimation with a neuralnetwork, 2017.
[31] A. Griffin. Conserving and gapless approximations for an inhomogeneous bose gas at finite
temperatures. Phys. Rev. B, 53:9341–9347, Apr 1996.
[32] T. L. Gustavson, P. Bouyer, and M. A. Kasevich. Precision rotation measurements with anatom interferometer gyroscope. Phys. Rev. Lett., 78:2046–2049, Mar 1997.
[33] M. N. Holger B. Nielsen. Dirac sea for bosons. Aug 1998.
[34] M. Holland, J. Williams, and J. Cooper. Bose-einstein condensation: Kinetic evolution obtained from simulated trajectories. Phys. Rev. A, 55:3670–3677, May 1997.
[35] M. J. Holland and K. Burnett. Interferometric detection of optical phase shifts at the heisenberg limit. Phys. Rev. Lett., 71:1355–1358, Aug 1993.
[36] J. Hu, L. Feng, Z. Zhang, and C. Chin. Quantum simulation of unruh radiation. Nature
Physics, 2019.
[37] N. M. Hugenholtz and D. Pines. Ground-state energy and excitation spectrum of a system
of interacting bosons. Phys. Rev., 116:489–506, Nov 1959.
[38] C. Jekeli. Navigation error analysis of atom interferometer inertial sensor. NAVIGATION,
52(1):1–14, 2005.
[39] C. Jekeli. Precision free-inertial navigation with gravity compensation by an onboard gradiometer. Journal of Guidance, Control, and Dynamics, 29(3):704–713, 2006.
[40] M. Kasevich and S. Chu. Atomic interferometry using stimulated raman transitions. Phys.
Rev. Lett., 67:181–184, Jul 1991.
[41] M. Khanahmadi and K. Mølmer. Time-dependent atomic magnetometry with a recurrent
neural network, 2020.
[42] D. Kingma and J. Ba. Adam: A method for stochastic optimization. International Conference
on Learning Representations, 12 2015.
[43] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing. Science,
220(4598):671–680, 1983.
[44] M. Kitagawa and M. Ueda. Squeezed spin states. Phys. Rev. A, 47:5138–5143, Jun 1993.
[45] S. J. J. M. F. Kokkelmans and M. J. Holland. Ramsey fringes in a bose-einstein condensate
between atoms and molecules. Phys. Rev. Lett., 89:180401, Oct 2002.
[46] S. J. J. M. F. Kokkelmans, J. N. Milstein, M. L. Chiofalo, R. Walser, and M. J. Holland.
Resonance superfluidity: Renormalization of resonance scattering theory. Phys. Rev. A,
65:053617, May 2002.

130[47] T. Kovachy, P. Asenbaum, C. Overstreet, C. A. Donnelly, S. M. Dickerson, A. Sugarbaker,
J. M. Hogan, and M. A. Kasevich. Quantum superposition at the half-metre scale. Nature,528(7583):530–533, Dec 2015.
[48] S. Kuang and S. Cong. Lyapunov control methods of closed quantum systems. Automatica,
44(1):98 – 108, 2008.
[49] Y. LeCun, Y. Bengio, and G. Hinton. Deep learning. nature, 521(7553):436, 2015.
[50] A. Lenef, T. D. Hammond, E. T. Smith, M. S. Chapman, R. A. Rubenstein, and D. E.
Pritchard. Rotation sensing with an atom interferometer. Phys. Rev. Lett., 78:760–763, Feb
1997.
[51] E. H. Lieb and W. Liniger. Exact analysis of an interacting bose gas. i. the general solution
and the ground state. Phys. Rev., 130:1605–1616, May 1963.
[52] L.-J. Lin. Self-improving reactive agents based on reinforcement learning, planning and teaching. Machine Learning, 8(3):293–321, May 1992.
[53] O. J. Luiten, M. W. Reynolds, and J. T. M. Walraven. Kinetic theory of the evaporative
cooling of a trapped gas. Phys. Rev. A, 53:381–389, Jan 1996.
[54] L. Mach. Ueber einen interferenzrefraktor, 1892.
[55] J. Mackeprang, D. B. R. Dasari, and J. Wrachtrup. A reinforcement learning approach forquantum state engineering. Quantum Machine Intelligence, 2(1):5, May 2020.
[56] P. J. Martin, B. G. Oldaker, A. H. Miklich, and D. E. Pritchard. Bragg scattering of atoms
from a standing light wave. Phys. Rev. Lett., 60:515–518, Feb 1988.
[57] A. A. Melnikov, H. Poulsen Nautrup, M. Krenn, V. Dunjko, M. Tiersch, A. Zeilinger, and
H. J. Briegel. Active learning machine learns to create new quantum experiments. Proceedingsof the National Academy of Sciences, 115(6):1221–1226, 2018.
[58] A. A. Melnikov, P. Sekatski, and N. Sangouard. Setting up experimental bell tests withreinforcement learning. Phys. Rev. Lett., 125:160401, Oct 2020.
[59] J. N. Milstein, C. Menotti, and M. J. Holland. Feshbach resonances and collapsing boseeinstein condensates. New Journal of Physics, 5:52–52, may 2003.
[60] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller. Playing atari with deep reinforcement learning, 2013.
[61] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves,
M. Riedmiller, A. K. Fidjeland, G. Ostrovski, et al. Human-level control through deep
reinforcement learning. Nature, 518(7540):529, 2015.
[62] H. P. Nautrup, N. Delfosse, V. Dunjko, H. J. Briegel, and N. Friis. Optimizing Quantum
Error Correction Codes with Reinforcement Learning. Quantum, 3:215, Dec. 2019.
[63] J. A. Nelder and R. Mead. A Simplex Method for Function Minimization. The Computer
Journal, 7(4):308–313, 01 1965.

131[64] J. Niederle, V. Popov, and L. Hlavatỳ. Functional Integrals in Quantum Field Theory andStatistical Physics. Mathematical Physics and Applied Mathematics. Springer Netherlands,
1983.
[65] M. Y. Niu, S. Boixo, V. N. Smelyanskiy, and H. Neven. Universal quantum control throughdeep reinforcement learning. npj Quantum Information, 5(1):33, Apr 2019.
[66] S. P. Nolan, A. Smerzi, and L. Pezzè. A machine learning approach to bayesian parameter
estimation, 2020.
[67] E. Oelker, R. B. Hutson, C. J. Kennedy, L. Sonderhouse, T. Bothwell, A. Goban, D. Kedar,
C. Sanner, J. M. Robinson, G. E. Marti, D. G. Matei, T. Legero, M. Giunta, R. Holzwarth,
F. Riehle, U. Sterr, and J. Ye. Demonstration of 4.8×10ˆ(-17) stability at 1 s for two independent optical clocks. Nature Photonics, 13(10):714–719, Oct 2019.
[68] M. Olshanii. Atomic scattering in the presence of an external confinement and a gas ofimpenetrable bosons. Phys. Rev. Lett., 81:938–941, Aug 1998.
[69] Z. Pagel, W. Zhong, R. H. Parker, C. T. Olund, N. Y. Yao, and H. Müller. Symmetric blochoscillations of matter waves. Phys. Rev. A, 102:053312, Nov 2020.
[70] R. H. Parker, C. Yu, W. Zhong, B. Estey, and H. Müller. Measurement of the fine-structure
constant as a test of the standard model. Science, 360(6385):191–195, 2018.
[71] L. H. Pedersen, N. M. Møller, and K. Mølmer. Fidelity of quantum operations. Physics
Letters A, 367(1):47 – 51, 2007.
[72] A. P. Peirce, M. A. Dahleh, and H. Rabitz. Optimal control of quantum-mechanical systems:
Existence, numerical approximation, and applications. Phys. Rev. A, 37:4950–4964, Jun 1988.
[73] M. E. Peskin and D. V. Schroeder. An Introduction to quantum field theory. AddisonWesley, Reading, USA, 1995.
[74] A. Peters, K. Y. Chung, and S. Chu. Measurement of gravitational acceleration by dropping
atoms. Nature, 400(6747):849–852, Aug 1999.
[75] D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov. Bose-einstein condensation in quasi-2d
trapped gases. Phys. Rev. Lett., 84:2551–2555, Mar 2000.
[76] H. F. Rice, V. Benischek, and L. Sczaniecki. Application of atom interferometric technology
for gps independent navigation and time solutions. In 2018 IEEE/ION Position, Location
and Navigation Symposium (PLANS), pages 1097–1106, 2018.
[77] F. Riehle, T. Kisters, A. Witte, J. Helmcke, and C. J. Bordé. Optical ramsey spectroscopy ina rotating frame: Sagnac effect in a matter-wave interferometer. Phys. Rev. Lett., 67:177–180,
Jul 1991.
[78] G. Rummery and M. Niranjan. On-line q-learning using connectionist systems. Technical
Report CUED/F-INFENG/TR 166, 11 1994.
[79] P. A. Ruprecht, M. J. Holland, K. Burnett, and M. Edwards. Time-dependent solution of
the nonlinear schrödinger equation for bose-condensed trapped neutral atoms. Phys. Rev. A,
51:4704–4711, Jun 1995.

132[80] L. Salvi, N. Poli, V. Vuletić, and G. M. Tino. Squeezing on momentum states for atominterferometry. Phys. Rev. Lett., 120:033601, Jan 2018.
[81] J. F. Schaff, T. Langen, and J. Schmiedmayer. Interferometry with atoms. 2015.
[82] C. Schubert, S. Abend, M. Gersemann, M. Gebbe, D. Schlippert, P. Berg, and E. M. Rasel.
Multi-loop atomic sagnac interferometry. Scientific Reports, 11(1):16121, Aug 2021.
[83] J. Schuff, L. J. Fiderer, and D. Braun. Improving the dynamics of quantum sensors withreinforcement learning. New Journal of Physics, 22(3):035001, mar 2020.
[84] A. Shankar, L. Salvi, M. L. Chiofalo, N. Poli, and M. J. Holland. Squeezed state metrology with bragg interferometers operating in a cavity. Quantum Science and Technology,
4(4):045010, oct 2019.
[85] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. van den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S. Dieleman, D. Grewe, J. Nham,
N. Kalchbrenner, I. Sutskever, T. Lillicrap, M. Leach, K. Kavukcuoglu, T. Graepel, and
D. Hassabis. Mastering the game of go with deep neural networks and tree search. Nature,529(7587):484–489, Jan 2016.
[86] D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre,
D. Kumaran, T. Graepel, T. Lillicrap, K. Simonyan, and D. Hassabis. A general reinforcement learning algorithm that masters chess, shogi, and go through self-play. Science,
362(6419):1140–1144, 2018.
[87] D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert,
L. Baker, M. Lai, A. Bolton, et al. Mastering the game of go without human knowledge.
Nature, 550(7676):354, 2017.
[88] K. Staliunas, S. Longhi, and G. J. de Valcárcel. Faraday patterns in bose-einstein condensates.
Phys. Rev. Lett., 89:210406, Nov 2002.
[89] R. Sutton, D. A. McAllester, S. Singh, and Y. Mansour. Policy gradient methods for reinforcement learning with function approximation. In NIPS, 1999.
[90] R. S. Sutton and A. G. Barto. Reinforcement learning: An introduction. MIT press, 2018.
[91] J. Taylor. Scattering theory: the quantum theory on nonrelativistic collisions. Wiley, 1972.
[92] G. Tesauro. Temporal difference learning and td-gammon. Commun. ACM, 38(3):58–68,
Mar. 1995.
[93] V. Vali and R. W. Shorthill. Fiber ring interferometer. Appl. Opt., 15(5):1099–1100, May
1976.
[94] H. Van Hasselt, A. Guez, and D. Silver. Deep reinforcement learning with double q-learning.
In Thirtieth AAAI conference on artificial intelligence, 2016.
[95] J. Wallnöfer, A. A. Melnikov, W. Dür, and H. J. Briegel. Machine learning for long-distancequantum communication. PRX Quantum, 1:010301, Sep 2020.
[96] D. F. Walls. Squeezed states of light. Nature, 306(5939):141–146, 1983.

133[97] R. Walser. Ground state correlations in a trapped quasi one-dimensional bose gas. Optics
Communications, 243(1):107 – 129, 2004. Ultra Cold Atoms and Degenerate Quantum Gases.
[98] Z. T. Wang, Y. Ashida, and M. Ueda. Deep reinforcement learning control of quantum
cartpoles. Phys. Rev. Lett., 125:100401, Sep 2020.
[99] C. J. C. H. Watkins and P. Dayan. Q-learning. In Machine Learning, pages 279–292, 1992.
[100] C. A. Weidner and D. Z. Anderson. Experimental demonstration of shaken-lattice interferometry. Phys. Rev. Lett., 120:263201, Jun 2018.
[101] C. A. Weidner and D. Z. Anderson. Simplified landscapes for optimization of shaken latticeinterferometry. New Journal of Physics, 20(7):075007, jul 2018.
[102] C. A. Weidner, H. Yu, R. Kosloff, and D. Z. Anderson. Atom interferometry using a shaken
optical lattice. Phys. Rev. A, 95:043624, Apr 2017.
[103] R. J. Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Mach. Learn., 8(3–4):229–256, May 1992.
[104] D. J. Wineland, J. J. Bollinger, W. M. Itano, F. L. Moore, and D. J. Heinzen. Spin squeezing
and reduced quantum noise in spectroscopy. Phys. Rev. A, 46:R6797–R6800, Dec 1992.
[105] L.-A. Wu, M. Xiao, and H. J. Kimble. Squeezed states of light from an optical parametricoscillator. J. Opt. Soc. Am. B, 4(10):1465–1475, Oct 1987.
[106] H. Xu, J. Li, L. Liu, Y. Wang, H. Yuan, and X. Wang. Generalizable control for quantum
parameter estimation through reinforcement learning. npj Quantum Information, 5(1):82,
Oct 2019.
[107] L. Zehnder. Ein neuer interferenzrefraktor, 1891.
[108] X.-M. Zhang, Z. Wei, R. Asad, X.-C. Yang, and X. Wang. When does reinforcement learning
stand out in quantum control? a comparative study on state preparation. npj Quantum
Information, 5(1):85, Oct 2019.
[109] H. Z. Zhigang Wu. Dynamics and density correlations in matter wave jet emission of a driven
condensate. Jun 2019.

Appendix A

Numerical methods for the Bose-Einstein condensate

A.1

Imaginary-time propagation
The imaginary-time propagation method is a numerical algorithm for finding the groundstate of the Schrödinger equation. It is particularly useful when the Schroödinger equation is nonlinear and exact diagonization of the Hamiltonian is difficult or not possible. This method involves
replacing the time t with the imaginary time −iτ in the time-dependent Schrödinger equation, i.e.,−∂Ψ= HΨ .
∂τThe solution to the equation above (we have set ℏ to 1 here) is Ψ(τ ) =(A.1)−Ei τ ψ , where Eiii ci ePand ψi are the eigenvalues and the eigenvectors respectively. At large τ , the state with the lowest
energy, i.e., the ground state, dominates the solution.
In practice, we start with an initial guess of the wavefunction, Ψ(0) =Pi ci ψi , and weneed to make sure that the ground state amplitude is non-zero, i.e., the guessed wavefunction is
not orthogonal to the ground state. We use numerical integration methods such as second order
Runge-Kutta (RK2) to advance in time, and the wavefunction needs to be renormalized every stepto avoid the loss of numerical precision. The total time it takes to reach the ground state is related
to the inverse of the energy gap between the ground state and the first excited state.
This method is analogous to the gradient descent method introduced in Chapter 3, with thecost function being the energy.

135

A.2

Spatial derivative in 1D
To calculate the spatial derivative in the kinetic term in the one dimensional Hamiltonian,

ℏ ∂
− 2m
, we use the finite differencing method, i.e.,
∂x2
2

2

ψ(x − ∆x) + ψ(x + ∆x) − 2ψ(x)
∂2ψ
≈
.
2
∂x
(∆x)2

(A.2)

The numerical error of this approach (central difference approximation) is O(∆x2 ). With theinfinite square well potential,

V (x) =





∞

if x ≤ 0 or x ≥ 1




0

if 0 < x < 1

,

(A.3)

we can impose the boundary condition, ψ(0) = 0 and ψ(1) = 0, and the wavefunction amplitudes
at these points stay zero throughout the simulation.

A.3

Spatial derivative in 2D
In Chapter 4, we considered a system in a disk geometry, and assumed that this system possess

cylindrical symmetry in two dimension, i.e., Ψ(r) = ψ(r). In this case, the spatial derivative in the
kinetic term in the Hamiltonian can be derived to be∇2 Ψ =∂ 2 ψ 1 ∂ψ
+
.
∂r2
r ∂r

(A.4)

We may use the finite differencing method outlined in the previous section to calculate the second
derivative. The central difference approximation for the first derivative is∂ψψ(r + ∆r) − ψ(r − ∆r)
≈
.
∂r
2∆r

(A.5)

We may also impose the boundary condition ψ(rmax ) = 0 at some cutoff radius. It is obvious thatat r = 0, Eq. (A.4) is ill-behaved with the factor 1/r, but physically one can expect that the kinetic
energy is still finite and remains a continuous function in space. Therefore, we need to find theappropriate boundary condition to mitigate the issue. The boundary condition at the center arises

136
from the symmetric assumption, ψ(r) = ψ(−r). This means that the first derivative is zero, andthat we can use the L’Hôpital’s rule for the second term in Eq. (A.4), which gives
1 ∂ψ
∂2ψ
=
.
r→0 r ∂r
∂r2
lim

(A.6)

With this boundary condition, the derivative in the kinetic term at the grid point r = 0 can becalculated as∇2 Ψ = 2

2 ψ(∆r) − 2 ψ(0)ψ(∆r) − ψ(0)
∂2ψ
≈2
=4
.
2
2
∂r
(∆r)
(∆r)2

(A.7)

Appendix B

Simulation techniques for the shaken lattice

B.1

Bloch’s theorem
Bloch’s theorem states that the wavefunction for each eigenstate in a spatially periodic systemwith periodicity a can be written as
ψ(x) = eiq·x uq (x) ,(B.1)

where q is the so-called quasi-momentum, and u(x) is a periodic function, i.e., u(x) = u(x + a).

B.2

Momentum basis
In the momentum basis, the Hamiltonian for the atoms in the optical lattice (defined in

Chapter 5 can be written as
Ĥ =

X p2
p



V0 iϕ
−iϕ
|p⟩ ⟨p| −e |p + 2ℏkL ⟩ ⟨p| + e|p⟩ ⟨p + 2ℏkL | .
2m
4

(B.2)

Note that the Hamiltonian only couples the momentum states with 2ℏkL spacing. As a result, themanifolds with different quasi-momentum q can be simulated independently. This is reminiscentof the Bloch’s theorem B.1 in the momentum representation. Assume that the initial population
distribution is peaked around the momenta p ∈ {2nℏkL |n ∈ Z}, and that the width of the peak issmaller than 2ℏkL , as can be expected from a wavefunction that spans across many lattice sites, thena convenient representation to use is |n, q⟩, which stands for the momentum state |p = 2nℏkL + ℏq⟩.
The narrow peak assumption implies that ⟨n′ , q ′ |n, q⟩ = δnn′ δ(q − q ′ ).

138
The trick mentioned above reduces the computation complexity from O(P 2 ) to O(N 2 Q),
where P represents the length of the momentum grid across the entire range that we are interested
in simulating, [2nmin + qmin , 2 nmax + qmax ] (factors of ℏ and kL are dropped), N is the length of thegrid for the momentum comb, N = nmax − nmin + 1, and Q is the length of the quasi-momentum
grid.

B.3

Accelerating reference frame
The Hamiltonian for a one-dimensional system in a reference frame accelerating at rate a in

the positive direction can be written as
Ĥ = Ĥ0 + max̂ ,(B.3)

where Ĥ0 is the original Hamiltonian in the rest frame. The second term acts as a fictitious potential
on an atom with mass m.
If the quantum state and Ĥ0 is represented in the momentum basis, as is the case in theprevious section, then we need to deal with the additional term with the position operator using
the split-step method. The approach we take is the symmetric Trotter decompositione(A+B)∆t ≈ eA∆t/2 eB∆t eA∆t/2 ,(B.4)

where the error is in fractions of ([A, [A, B] + 2[B, [A, B])∆t3 /24 . We associate A and B with theoperators −iĤ0 /ℏ and −imax̂/ℏ respectively. The evolution e−iĤ0 ∆t/ℏ can be simulated using theRK2 method with the formulation we set up in the previous section. The evolution that involves
the position operator can be simulated as a shift operator, i.e.,e−imax̂∆t/ℏ |p⟩ = |p − ma∆t⟩ .

(B.5)

Since this shift operator acts the same on all |n, q⟩ states, we can simply shift the whole quasimomentum grid, q → q − ma∆t/ℏ. This trick avoids the need to perform fast Fourier transforms
when switching between the momentum and position evolutions, which are computationally expensive.

ProQuest Number: 29258937
INFORMATION TO ALL USERS
The quality and completeness of this reproduction is dependent on the quality
and completeness of the copy made available to ProQuest.

Distributed by ProQuest LLC ( 2022 ).
Copyright of the Dissertation is held by the Author unless otherwise noted.

This work may be used in accordance with the terms of the Creative Commons licenseor other rights statement, as indicated in the copyright statement or in the metadata
associated with this work. Unless otherwise specified in the copyright statement
or the metadata, all rights are reserved by the copyright holder.

This work is protected against unauthorized copying under Title 17,
United States Code and other applicable copyright laws.

Microform Edition where available © ProQuest LLC. No reproduction or digitizationof the Microform Edition is authorized without permission of ProQuest LLC.

ProQuest LLC
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, MI 48106 - 1346 USA