University of Cambridge
Department of Applied Mathematics
and Theoretical Physics (DAMTP)
and
Wolfson College
Thermodynamic and hydrodynamic
behaviour of interacting Fermi gases
Olga Goulko
Dissertation submitted for the degree of Doctor of Philosophy
November 10, 2011

Declaration
This dissertation is the result of my own work and includes nothing which is the outcome
of work done in collaboration, except where specifically indicated in the text. No part of
this dissertation has been previously submitted for a degree or any other qualification.
Olga Goulko
November 10, 2011
2
Summary
Fermionic matter is ubiquitous in nature, from the electrons in metals and semiconduc-
tors or the neutrons in the inner crust of neutron stars, to gases of fermionic atoms, like
40K or 6Li that can be created and studied under laboratory conditions. It is especially
interesting to study these systems at very low temperatures, where we enter the world
of quantum mechanical phenomena. Due to the Fermi-Dirac statistics, a dilute system
of spin-polarised fermions exhibits no interactions and can be viewed as an ideal Fermi
gas. However, interactions play a crucial role for fermions of several spin species.
This thesis addresses several questions concerning interacting Fermi gases, in par-
ticular the transition between the normal and the superfluid phase and dynamical
properties at higher temperatures. First we will look at the unitary Fermi gas: a
two-component system of fermions interacting with divergent scattering length. This
system is particularly interesting as it exhibits universal behaviour. Due to the strong
interactions perturbation theory is inapplicable and no exact theoretical description is
available. I will describe the Determinant Diagrammatic Monte Carlo algorithm with
which the unitary Fermi gas can be studied from first principles. This algorithm fails in
the presence of a spin imbalance (unequal number of particles in the two components)
due to a sign problem. I will show how to apply reweighting techniques to generalise
the algorithm to the imbalanced case, and present results for the critical temperature
and other thermodynamic observables at the critical point, namely the chemical po-
tential, the energy per particle and the contact density. These are the first numerical
results for the imbalanced unitary Fermi gas at finite temperature. I will also show
how temperatures beyond the critical point can be accessed and present results for the
equation of state and the temperature dependence of the contact density.
At sufficiently high temperatures a semiclassical description captures all relevant
physical features of the system. The dynamics of an interacting Fermi gas can then
be studied via a numerical simulation of the Boltzmann equation. I will describe such
a numerical setup and apply it to study the collision of two spin-polarised fermionic
clouds. When the two components are separated in an elongated harmonic trap and
then released, they collide and for sufficiently strong interactions can bounce off each
other several times. I will discuss the different types of the qualitative behaviour, show
how they can be interpreted in terms of the equilibrium properties of the system, and
explain how they relate to the coupling between different excitation modes. I will also
demonstrate how transport coefficients, for instance the spin drag, can be extracted
from the numerical data.
3
Acknowledgements
First and foremost thanks go to my Ph.D. supervisor Matthew Wingate who has intro-
duced me to this interesting topic and has guided me through my research. Without
his help, advice and encouragement this thesis would not have been possible. Matt has
always found time to answer all my questions and to explain things to me that I did not
understand. I could not have hoped for a more approachable and patient supervisor.
I would also like to thank my other collaborators, Fre´de´ric Chevy and Carlos Lobo.
I very much enjoyed our inspiring discussions and the great working atmosphere. I hope
that in the future we will continue to work together.
Throughout my Ph.D. many people from the department have been a great help.
I am especially grateful to Francis Bursa, Johannes Hofmann, Ron Horgan, Adrian
Kent, Nick Manton and Stefan Meinel for many interesting discussions and helpful
advice. Thanks go to all members of the Lattice Field Theory group for the friendly
atmosphere and the interesting seminars and discussions. A big thank you also to
Amanda Stagg for help with numerous organisational issues.
Many colleagues at other institutions have been a tremendous help and provided
me with useful advice and new ideas. I am very grateful to Evgeni Burovski, Yvan
Castin, Joaqu´ın Drut, Kris van Houcke, Dean Lee, Meera Parish, Nikolay Prokof’ev,
Christophe Salomon, Ariel Sommer, Boris Svistunov, Michael Urban, Herbert Wagner,
Felix Werner and Martin Zwierlein. Special thanks go to Evgeni Burovski and Felix
Werner for providing me with their numerical data which has been included in this
thesis. Numerous other people whom I met during conferences and seminars have given
me new ideas and helped me to learn new things. In particular I would like to thank
all participants of the ECT* workshop “Sign problems and complex actions” for many
interesting talks and discussions.
I cannot express enough thanks to my family and friends for their constant encour-
agement and moral support during my time in Cambridge. A huge thank you goes to
my fellow Part III and Ph.D. students at DAMTP for many interesting and motivating
discussions and study group meetings. Without the wonderful friendly atmosphere at
Wolfson College, the inspiring lunch conversations and the occasional thesis “support
group” meetings this work could not have been completed.
My research was funded by the German Academic Exchange Service (DAAD), the
Engineering and Physical Sciences Research Council (EPSRC) and the Cambridge Eu-
ropean Trust, to all of whom I am very grateful for their generous financial support.
I would also like to thank the Department of Applied Mathematics and Theoretical
4
Physics (DAMTP), Wolfson College and the Cambridge Philosophical society for pro-
viding me with travel grants that have enabled me to attend numerous workshops and
conferences. For local financial support during conferences I am grateful to the ECT*
in Trento, Durham University and the ENS in Paris. Many thanks go to all the seminar
organisers who have invited me to give talks at their institutions. I would like to thank
the Institute for Nuclear Theory at the University of Washington for its hospitality and
financial support during the programme “Fermions from Cold Atoms to Neutron Stars:
Benchmarking the Many-Body Problem”.
This research was very computationally intensive and would not have been possible
without the support and the helpful advice of the computing officers at DAMTP, es-
pecially Elie Bassouls, Mike Rose and John Sutton. I am also very grateful to Stefan
Meinel and Hanno Rein for their help with various programming issues. This work has
made use of the resources provided by the University of Cambridge High Performance
Computing Service (http://www.hpc.cam.ac.uk/). I am grateful to Stuart Rankin
from the HPC for help on various occasions.
5
Contents
Summary 3
Acknowledgements 4
1 Introduction 8
1.1 Ultracold atomic gases . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2 Interacting Fermi gases . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 The unitary Fermi gas . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.4 Spin imbalance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.5 Dynamical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2 The Fermi-Hubbard model at finite temperature 22
2.1 Dispersion relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2 Tuning the coupling constant . . . . . . . . . . . . . . . . . . . . . . . 24
2.3 Lattice corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4 The finite temperature formalism . . . . . . . . . . . . . . . . . . . . . 29
2.5 Matrix notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3 Order parameter and finite-size scaling 33
4 Diagrammatic Monte Carlo sampling 39
4.1 Formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2 Extending the configuration space . . . . . . . . . . . . . . . . . . . . . 41
4.3 Sampling with a different distribution . . . . . . . . . . . . . . . . . . . 42
5 Implementation of the algorithm 44
6 Observables at the critical point 50
6.1 Spin susceptibility and pseudogap . . . . . . . . . . . . . . . . . . . . . 51
6.2 Critical temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
6.3 Energy per particle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6.4 Chemical potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.5 Contact density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.6 Comparison with the literature . . . . . . . . . . . . . . . . . . . . . . 63
6
7 Observables beyond the critical point 67
7.1 The equation of state for the density . . . . . . . . . . . . . . . . . . . 69
7.2 The equation of state for the pressure . . . . . . . . . . . . . . . . . . . 71
7.3 The temperature dependence of the contact . . . . . . . . . . . . . . . 74
8 Simulation of the Boltzmann equation 77
8.1 Distributions in a trap . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
8.2 The Boltzmann equation . . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.3 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
8.4 Tests and optimal parameters . . . . . . . . . . . . . . . . . . . . . . . 87
9 Collision of two fermionic clouds 95
9.1 Qualitative Behaviour . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
9.2 Transitions between regimes . . . . . . . . . . . . . . . . . . . . . . . . 97
9.3 Coupling between excitation modes . . . . . . . . . . . . . . . . . . . . 98
9.4 Transport coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
10 Conclusions and Outlook 109
A Matrix update procedure 112
A.1 The Sherman-Morrison Formula . . . . . . . . . . . . . . . . . . . . . . 112
A.2 Adding one row and one column . . . . . . . . . . . . . . . . . . . . . . 114
A.3 Removing one row and one column . . . . . . . . . . . . . . . . . . . . 117
A.4 Changing the elements of one row . . . . . . . . . . . . . . . . . . . . . 123
Bibliography 126
7
Chapter 1
Introduction
1.1 Ultracold atomic gases
The first experimental realisations of a Bose-Einstein condensate (BEC) in a dilute
atomic gas in 1995 marked the birth of what has by now become an entire new area
of research within condensed matter physics [1, 2]. Since then enormous progress has
been made both theoretically and experimentally in the field of cold atoms. The typical
densities of the systems considered range between 1013 and 1015 cm−3, which is four to
six orders of magnitude lower than of the air surrounding us, not to mention denser
systems like liquids or solids. To observe quantum phenomena at such low densities
temperatures below 10−5 K are required – a major technical challenge. Only with the
development of advanced cooling techniques like laser cooling and evaporative cooling
could such low temperatures be achieved.
But why is the study of ultracold gases so important in the first place? Partial
Bose-Einstein condensation can be observed in liquid 4He already at temperatures be-
low 2.17 K, which is why over five decades passed between the discovery of superfluidity
in liquid helium and the first realisation of a gaseous BEC. The crucial difference be-
tween the two systems is that 4He is a liquid, which means that strong correlations are
induced by interatomic interactions. These complicated interactions obscure a mean-
field description for conventional superfluids and superconductors. Dilute gases on the
other hand can form essentially pure condensates, well-described by the macroscopic
wave function. Hence they can provide a direct realisation of many basic models in
condensed matter physics.
The significance of cold atomic gases can be summarised in two main points: good
tunability and a multitude of analogues. Dilute gases can be easily manipulated with
lasers and magnetic fields, which gives experimentalists excellent control over a vari-
ety of properties, such as interaction strength, temperature, trapping potential and
even dimension. This remarkable tunability allows to enter regimes that have not been
experimentally accessible before. Also, thanks to the low density of the system, micro-
scopic length scales are large enough to be easily observed by optical means. Superfluid
features, collective oscillations, quantised vortices, few-particle systems in microtraps
and Mott-insulator to superfluid transitions in optical lattices are just a small selection
8
of the many phenomena that have so far been subject of research in this field [3, 4].
At low temperatures the thermal de Broglie wavelength is large compared to the
range of the interatomic potential and the exact form of the interaction becomes unim-
portant for the macroscopic properties of the system. Therefore the model of a dilute
quantum gas can represent various systems in nature and can be used in different areas
of physics. For example it finds application in nuclear physics where results from the
field of cold gases can be used to describe systems such as low-density neutron mat-
ter or even atomic nuclei. The neutrons in the core of a neutron star are believed to
form a fermionic condensate that can be studied with the methods of cold Fermi gases.
Even quantum chromodynamics (QCD) at high temperatures and densities might have
similar physics to these models.
In this thesis I will focus on different features of interacting Fermi gases, ranging
from superfluidity to spin transport. I will begin with a general overview on interacting
Fermi gases introducing the main theoretical concepts. Chapters 2-6 are devoted to a
diagrammatic Monte Carlo study of the two-component unitary Fermi gas at the critical
point. This work was done in collaboration with Matthew Wingate and results were
published in [5], [6] and [7]. Chapter 7 extends this study to temperatures beyond the
critical point, which is still work in progress. In the chapters 8 and 9 I will describe how
different dynamical properties of a trapped interacting Fermi gas at high temperatures
and finite values of the scattering length can be studied with a numerical simulation
of the Boltzmann equation. This work was done in collaboration with Fre´de´ric Chevy
and Carlos Lobo and some of the results appeared in [8]. A summary of all results and
an outlook on future work can be found in chapter 10. Throughout the text we will use
units in which ~ = kB = 1.
1.2 Interacting Fermi gases
The physical consequences of quantum degeneracy are very different for Bose and Fermi
gases [3, 4]. Unlike an ideal Bose gas, which exhibits a phase transition to a BEC at low
temperatures, for a non-interacting Fermi gas quantum degeneracy only corresponds to
a smooth crossover from classical to quantum behaviour. This is due to Pauli’s exclu-
sion principle that does not allow two identical fermions to occupy the same quantum
state and thus inhibits s-wave scattering for non-interacting fermions. The existence
of a phase transition can only be due to interactions. Depending on the depth of the
interaction potential the system exhibits a rich variety of features.
In dilute systems interparticle interactions can be accurately described through two-
body collisions [1, 9, 10]. Assuming elastic collisions and neglecting small relativistic
spin interactions the problem of a collision between two particles of equal mass m
reduces to the solution of the Schro¨dinger equation for the scattering of one particle
9
with reduced mass m/2 and positive energy E = k2/m > 0 on the spherically symmetric
finite range potential V (r),
(
−∇
2
m + V (r)− E
)
ψ(r) = 0. (1.1)
At large distances from the scatterer (r À r0, where r0 is the range of the potential) the
solution of the Schro¨dinger equation is a superposition of the original incoming plane
wave propagating in the direction of k with an outgoing spherical wave,
ψ(r) = eik·r + f(k, θ)e
ikr
r , (1.2)
where the scattering amplitude f(k, θ) depends on the absolute value of the wavevector
k and the angle θ between the directions of r and k. The scattering amplitude tells
us about the probability for the scattered particle to pass through the surface element
r2dΩ per unit time, ∣∣∣∣f
eikr
r
∣∣∣∣
2
vr2dΩ = v|f |2dΩ, (1.3)
where dΩ is the solid angle element. The quantity dσ = |f |2dΩ has the dimension
of an area and is called the differential cross section. An explicit expression for the
scattering amplitude can be obtained by solving the Schro¨dinger equation (1.1) for
r À r0 and matching it to the asymptotic behaviour (1.2). The spherical symmetry of
the potential implies that the solution must be axially symmetric around the direction
k of the incoming flux. Hence we can decompose the wave function into the `-wave
components of angular momentum,
ψ(r) =
∞∑
`=0
P`(cos θ)χk`(r)kr , (1.4)
where P` are the Legendre polynomials and the functions χk`(r) satisfy the radial
Schro¨dinger equation,
d2χk`
dr2 −
`(`+ 1)
r2 χk` +m(E − V (r))χk` = 0. (1.5)
To match to the asymptotic form (1.2) we note that the plane wave eik·r also has a
spherical wave expansion [10],
eik·r = eikr cos θ = 12ikr
∞∑
`=0
(2`+ 1)P`(cos θ)
(eikr − e−i(kr−pi`)) . (1.6)
As the Legendre polynomials are pairwise orthogonal we need to match the different
` contributions individually. In the low-energy limit of ultracold collisions (k ¿ 1/r0)
10
the only relevant contribution comes from the s-wave component ` = 0. In this case
the radial Schro¨dinger equation outside the range of the potential r À r0 simplifies to
d2χk0
dr2 +mEχk0 = 0 (1.7)
and has the solution
χk0 = A0 sin(kr + δ0), (1.8)
with two integration constants: A0 and the phase shift δ0, which generally depend on
k. Comparing with (1.2) we obtain A0 = eiδ0 and
f`=0(k) = e
2iδ0 − 1
2ik =
1
k cot δ0 − ik . (1.9)
The ` = 0 contribution to the scattering amplitude turns out to be independent of the
angle θ, hence s-wave scattering is equally likely in any direction. Taking the limit
k → 0 we can neglect the energy term altogether, and the solution of (1.7) simplifies to
a linear function,
χk0(r) ∝ r − a, (1.10)
where the coefficient a has the dimension of length and is referred to as the scattering
length. The value (and sign) of the scattering length depends on the potential V (r) and
can be calculated by solving the Schro¨dinger equation inside the range of the potential
and matching it to the asymptotic form (1.10) at r = r0. The scattering length is
then the intercept of this asymptote. As an example Fig. 1.1 shows the solution of the
Schro¨dinger equation at small positive energies for the square potential
V (r) =
{
V0 if r ≤ r0
0 if r > r0
for different values of V0. Inside the range of this potential the solution is a sine function
if V0 < 0 (potential well) and an exponential function if V0 > 0 (potential barrier). We
can see from Fig. 1.1(a) that the scattering length is positive for a repulsive potential
with V0 > 0 and from Fig. 1.1(b) that it is negative for a weakly attractive potential
with V0 < 0. As we increase the potential depth the scattering length approaches −∞
and for even deeper attractive potentials becomes positive again, as shown in Fig. 1.1(c).
By comparing (1.10) to the k → 0 limit of (1.8) we can obtain the leading order
expression for the phase shift δ0,
lim
k→0
χk0 ∝ limk→0 sin
(
k
(
r + δ0k
))
= r + δ0k (1.11)
11
V0
0 r0 ra>0
(a) repulsive
V0
0 r0 ra<0
(b) weakly attractive
V0
0 r0 ra>0
(c) strongly attractive
Figure 1.1: The radial wave function and its asymptote for different values of potential depth
V0. The scattering length a is the intercept of the asymptote with the r-axis.
and therefore the scattering length equals a = −δ0/k. Hence the phase shift δ0 is
proportional to k. Higher ` phase shifts only contribute to order k2`+1 and are thus
indeed negligible at low energies [9, 10]. For the s-wave scattering amplitude (1.9) this
implies
f`=0(k) = − 1a−1 + ik , (1.12)
to leading order in k. Note that when k → 0 the s-wave scattering amplitude tends
to a constant value, f`=0(k → 0) = −a. Hence the scattering length alone is sufficient
to describe low-energy scattering. In particular, the exact shape of the interaction
potential becomes unimportant. If needed, we can revert to an effective potential
description and replace the potential by a more convenient form, for instance a zero-
range potential, provided that it has the same long-distance behaviour. For the s-wave
cross section we obtain
σ =
∫
dσ =
∫
|f |2dΩ =
∫ sin θdθdφ
a−2 + k2 =
4pia2
1 + k2a2 . (1.13)
Note that for identical fermions the orbital part of the wave function must be anti-
symmetric and hence the cross section vanishes, dσ = |f(θ) − f(pi − θ)|2dΩ = 0. This
explains why s-wave scattering is inhibited for identical fermions and a spin polarised
system behaves like an ideal gas. To observe low-energy interactions at least two dif-
ferent fermionic components are required.
Including the next-to-leading order in the expansion of the phase shift we obtain
k cot δ0 = −a−1 + 12rek
2, (1.14)
where re is called the effective range of the potential. Due to time reversal symmetry
no linear term in k is present [11]. Note that the effective range is generally not the
same as the range of the potential r0 and, like the scattering length, can have arbitrary
sign.
12
∞ −∞0
1
a
BEC regime
strongly bound (bosonic)
molecules of two fermions
Unitarity
strongly interacting
fermions
BCS regime
pairs of fermions weakly
bound in momentum space
Figure 1.2: The different regimes of behaviour as a function of the inverse scattering length
So far we have only considered positive energy solutions. If the energy is negative,
solutions of the Schro¨dinger equation are two-body bound states. The asymptotic
(r À r0) bound state wave function for ` = 0 has the form
χb(r) ∝ e−
√
m|Eb|r, (1.15)
where Eb < 0 is the binding energy of the dimer and both particles are again assumed
to have equal mass m. Expanding this expression for small energies Eb to leading order,
χb(r) ∝ 1−
√
m|Eb|r, (1.16)
and comparing the result to equation (1.10) yields a−1 =
√
m|Eb| > 0 and hence
Eb = −1/ma2. So bound states occur at positive scattering length, with binding energy
inversely proportional to the square of the scattering length. This implies that a small
positive scattering length can correspond both to a repulsive potential, as depicted in
Fig. 1.1(a) and a bound state with large binding energy. As we shall see, the repulsive
Fermi gas is unstable, and it is the gas of dimers that is usually realised experimentally.
To sum up, we have seen that in the presence of interactions the behaviour of a
system depends crucially on sign and value of the scattering length a. In the following
we will always consider fermions of two species, for instance two different spin states, as
otherwise no s-wave interactions are possible. As sketched in Fig. 1.2, we can classify
the interactions into three qualitative regions: small and negative a (the BCS regime),
small and positive a (the BEC regime) and the unitary regime of divergent a that, in
terms of the inverse scattering length, lies in between the BCS and the BEC regime.
The regime of small and negative scattering length, a < 0 and kF |a| ¿ 1, (where
13
kF = (3pi2n)1/3 is the Fermi wavevector defined in terms of the total number density n)
corresponds to a weakly attractive gas and can be accurately described using the theory
J. Bardeen, L. N. Cooper and J. R. Schrieffer introduced to explain superconductivity
[3]. It is thus referred to as the BCS regime. The BCS theory assumes the existence of
a weakly attractive force between the fermions. Below some critical temperature, this
force will cause atoms with opposite momenta on the Fermi surface to form Cooper
pairs – bound states in momentum space with zero centre of mass momentum and
exponentially small binding energy. This mechanism is similar to the formation of
Cooper pairs of electrons in a superconducting metal. These weakly bound pairs are
very large and overlap strongly in coordinate space. They form a condensate in the zero-
momentum state, which gives rise to superfluidity, analogous to the superconducting
state of certain metals. This regime is well-understood and the BCS theory provides
exact results for many interesting quantities, both at zero and finite temperature. The
critical temperature for example is
Tc ≈ 0.28TF epi/2kF a, (1.17)
where TF = εF = k2F/2m stands for the Fermi temperature. With decreasing kF |a| the
critical temperature becomes exponentially small, which makes the phase transition
difficult to study experimentally.
As we have seen above, bound states of two atoms can be formed at positive scatter-
ing length a, with energy Eb = −1/ma2 that becomes large for small values of a. Such
strongly bound dimers are bosons since they consist of two fermions. They behave like
a Bose gas and form a BEC below a critical temperature, so that the regime a > 0,
kF |a| ¿ 1 is referred to as the BEC regime. A theoretical description of the BEC
regime is provided by the usual theory for Bose gases. The critical temperature is given
by
Tc = 0.218TF , (1.18)
and is considerably higher than in the BCS case, which makes it much easier to observe
in an experiment. It is important to stress that the bound states in the BEC regime are
strongly bound molecules, and not weakly bound pairs as in the BCS case. These two
regimes are fundamentally different and can be distinguished in an experiment, since in
the BEC regime the formation of molecules is followed by their condensation only when
the molecule density becomes high enough, while in the BCS regime pair formation and
condensation take place simultaneously.
With respect to the relevant degrees of freedom – fermionic pairs in the BCS case
and bosonic molecules in the BEC regime – these two limits correspond to weakly
interacting systems and can be described analytically. The region in between, the
so-called BEC-BCS crossover, is more challenging. In the limit a → ∞, the s-wave
14
scattering amplitude obeys the universal law f`=0(k) = i/k, independent of the form
of the interaction. This is referred to as the unitary Fermi gas, which we will discuss
in detail in section 1.3. In this case f`=0(k) takes its maximal value – the system is
maximally interacting.
Experimentally, the absence of s-wave scattering in a gas of identical fermions leads
to major difficulties with evaporative cooling, a technique based on scattering. This
difficulty can be bypassed with help of “sympathetic cooling” techniques, which either
employ two different spin components of the same Fermi gas or add a Bose gas com-
ponent to the Fermi gas as a refrigerant. Another important tool for the experimental
study of Fermi gases interacting with different values of the scattering length is provided
by Feshbach resonances, which allow to vary the value and even the sign of the scat-
tering length by simply tuning a magnetic field [3, 4, 12]. A Feshbach resonance occurs
when the energy associated with the scattering between two particles (open channel)
becomes very close to the energy of a bound state relative to a different channel (closed
channel). The scattered particles can then temporarily enter the metastable bound
state, which can significantly increase the scattering cross section. If open and closed
channel carry different magnetic momenta, their energies can be continuously shifted
with respect to each other via the Zeeman effect in the hyperfine levels, by varying
an external uniform magnetic field. This can be used to tune the energies towards
the resonance. By changing the energy of the closed channel from just below to just
above the energy in the open channel in the presence of a small coupling, the scattering
length can be tuned from large and positive to large and negative values. Crossing this
threshold, the scattering length diverges. The dependence of the scattering length on
the external magnetic field B can be phenomenologically modelled by the relation
a = abg
(
1− ∆BB −B0
)
, (1.19)
where abg is the background scattering length for collisions entirely in the open channel,
and ∆B and B0 describe the width and position of the resonance respectively [4, 11].
Such a resonance for the two lowest magnetic sub-states of 6Li is plotted in Fig. 1.3.
The whole range of possible values for the scattering length a, including unitarity, can
be experimentally achieved using Feshbach resonances.
As we have seen, small and positive values of the scattering length can correspond to
both a repulsive gas and a gas of strongly bound dimers. When the scattering length is
tuned starting from a Feshbach resonance the atoms are captured in bound states, and
hence only the latter case is realised [13]. It is also possible to experimentally create a
repulsive gas by adiabatically switching on the scattering length starting from a = 0.
The repulsive Fermi gas is unstable and does not exhibit superfluidity.
15
-20
-15
-10
-5
 0
 5
 10
 15
 20
 600  700  800  900  1000  1100
sc
att
eri
ng
 le
ng
th 
(10
00a
0)
magnetic field (G)
Figure 1.3: Magnetic field dependence of the scattering length (1.19) for the two lowest
magnetic sub-states of 6Li. The Feshbach resonance (dashed line) occurs at
B0 = 834 G with a width ∆B = 300 G. The background scattering length
abg = −1405a0 is very large, a0 is the Bohr radius.
1.3 The unitary Fermi gas
The regime of divergent scattering length a →∞ is called unitarity and is particularly
interesting for a number of reasons. In this case the gas is both dilute (the range of
the interatomic potential is much smaller than the interparticle distance) and strongly
interacting (the scattering length is much larger than the interparticle distance) at the
same time [3]. As discussed in the previous section, details of the interatomic potential
are unimportant for the description of low-energy interactions and the scattering length
is the only relevant interaction-related scale of the problem. At unitarity the scattering
length diverges so that this length scale is no longer present. The behaviour of the gas
now only depends on two dimensionful parameters: the temperature and the density of
the system. Hence all thermodynamic observables must be universal functions of the
temperature T and the Fermi energy of the non-interacting gas εF = (3pi2n)2/3/2m.
For instance the critical temperature is simply a dimensionless number times the Fermi
energy, while the chemical potential in units of the Fermi energy is a purely universal
function of the temperature.
On the practical side the unitary Fermi gas offers fascinating prospects to study high-
temperature superfluidity. Compared to other known systems the critical temperature
of the phase transition into the superfluid phase is remarkably high, as depicted in
Table 1.1. If we scale it to the density of electrons in a metal this transition would occur
far above room temperature [11]. Exploring the properties of the unitary Fermi gas
promises valuable insights into the mechanisms behind high-temperature superfluidity
and superconductivity. Another advantage of the unusually high critical temperature
is that the phase transition becomes experimentally accessible, as opposed to the BCS
16
Tc TF Tc/TF
metallic superconductors 1 – 10 K 50 – 150 kK 10−5 – 10−4
superfluid 3He 2.6 mK 5 K 5 · 10−4
high-temperature superconductors 35 – 140 K 2000 – 5000 K 10−2
neutron stars 1010 K 1011 K 10−1
strongly interacting atomic Fermi gases 200 nK 1 µK 0.2
Table 1.1: Typical critical temperatures Tc, Fermi temperatures TF , and their ratios for
various fermionic superfluids and superconductors, see [11].
regime where Tc is exponentially small.
Due to the strength of the interaction and the consequent lack of a small param-
eter for perturbation theory, no exact theory for the Fermi gas at large values of the
scattering length has been developed yet. There exist several different analytical and
numerical approaches. One possibility is to generalise the BCS mean-field theory to all
values of the scattering length [14]. However at finite temperature fluctuations become
crucial and the mean-field approximation is no longer valid. And even at zero tempera-
ture the many-body problem along the BCS-BEC crossover cannot be solved, since the
BCS model only contains zero-momentum pairs, and thus density fluctuations are not
accounted for [4].
One example of an analytical approach is to postulate N fermionic flavours and to
consider a system with 2N (rather than only two) fermionic fields, interacting with a
bosonic field [15]. The fermionic degrees of freedom can be integrated out, leaving only
the complex scalar field with a corresponding effective action. From this an expan-
sion in 1/N can be generated for different observables. The problem is solved in the
N → ∞ limit with correction terms of leading order in 1/N . The 1/N corrections are
comparatively large, given that the physically interesting case is N = 1.
A self-consistent many-body approach based on the Luttinger – Ward formalism is
presented in [16]. In this approach the grand thermodynamic potential is expressed
in terms of Green’s functions and interactions between fermions are represented as a
perturbative series of irreducible Feynman diagrams. As the complete series is impos-
sible to calculate one needs to employ the ladder approximation which describes the
formation of pairs and is exact in the BCS limit. The De Dominicis – Martin extention
allows to also incorporate bosonic degrees of freedom and improves the results across
the BEC-BCS crossover. An advantage of this approach is that the Green’s function
need not be approximated in terms of free Green’s functions. Also the results are valid
for all values of the scattering length. However it fails to correctly capture the critical
behaviour at the superfluid transition, as it predicts a weak first order transition rather
than a second order one.
Other analytical methods include an effective theory based on an expansion around
d = 4 − ² and d = 2 − ² spatial dimensions [17, 18]. Both limits can be treated ana-
17
lytically: near d = 4 the unitary Fermi gas behaves like a weakly interacting system
of fermionic and bosonic degrees of freedom and near d = 2 like a weakly interacting
Fermi gas, see [19] for a review. The physical case d = 3 is obtained through extrap-
olation of the series expansion to ² = 1. To first order, this approach allows to derive
an effective potential, from which thermodynamical quantities can be calculated with
good accuracy. However, the next order reveals a deviation from the numerically and
experimentally predicted values, which indicates that the expansion might not converge
sufficiently quickly. A calculation to third or higher order involves a significant increase
in difficulty.
Due to the strong nonperturbative nature of the problem numerical approaches are
the only ones that can give reliable quantitative predictions about the properties of the
Fermi gas in the unitarity limit. Unlike the analytical approximations they can start
directly from first principles and model the system in a systematically improvable way.
They provide quantitative results that can be compared directly with experimental
data and can be used as a benchmark to test analytical methods. An example is the
beyond mean-field approximation from [20], where the equation of state of the unitary
Fermi gas is compared to Monte Carlo results using three different perturbation (T -
matrix) schemes. Since there is no a priori judgement as to which perturbation scheme
is more accurate, the comparison to numerical results is vital. Monte Carlo results are
also used for a phenomenological interpolation of the analytical prediction around the
critical point, where the purely analytical results are poor.
By today, numerical methods have advanced so far that they can produce high
precision results for several quantities, for instance the critical temperature, with up to
a few percent accuracy. Many numerical approaches use methods from nuclear theory,
for instance the effective field theory approach in [21] which uses hybrid Monte Carlo
to simulate cold dilute neutron matter on the lattice. At zero temperature fixed-node
diffusion Monte Carlo finds wide application. In this setup the wave function is sampled
via a trial function using the Schro¨dinger equation. The fixed node approximation
enforces the positive definiteness of the sampling function, to avoid the fermionic sign
problem. At finite temperature it is used for the nonperturbative restricted path integral
Monte Carlo [22]. Due to the fixed-node approximation at finite temperature this
approach can only give upper bounds on the critical temperature. The auxiliary field
quantum Monte Carlo method [23, 24, 25] uses a continuum model with a momentum
cut-off and is independent of the dispersion relation.
In this work we will focus on the determinant diagrammatic Monte Carlo (DDMC)
approach [26], which will be described in detail in the following chapters. This approach
avoids the fermionic sign problem for the two-component balanced gas and can provide
precise predictions for the critical temperature and other thermodynamic observables
at unitarity.
18
1.4 Spin imbalance
So far, most numerical studies were limited to the balanced case, when the number
of fermions in the two spin components is equal. An imbalance in the particle num-
ber results in new interesting effects, which makes a detailed numerical study desirable
[3, 27, 28]. Amongst the main areas of interest are the finite temperature phase dia-
gram, new exotic phases and phase transitions, as well as the analysis of the strongly
interacting normal state at high imbalance. Imbalanced Fermi gases are also related
to superconductors under an external magnetic field. Since the orbital motion of the
electrons screens this field in bulk superconductors, a highly imbalanced state is difficult
to achieve. Superfluid Fermi gases do not suffer from this restriction. Also, unlike with
superconductors, we are not limited to the BCS regime but can work across the whole
BEC-BCS crossover.
As the single component Fermi gas does not interact via s-wave scattering and
consequently cannot undergo a phase transition into the superfluid phase it is clear that
as we increase the imbalance at some point superfluidity must break down even at zero
temperature. Experiments agree however that the superfluid phase is remarkably stable
towards imbalance and that this breakdown occurs at a very high critical imbalance of
around 77% at unitarity [29, 30]. The authors of [31] quote an even higher value of over
90%, but this discrepancy is probably due to their system being out of equilibrium, as a
consequence of the highly elongated trap geometry and the evaporative cooling scheme
which create a polarisation pattern that favours a superfluid at the trap centre [32]. In
all experiments a phase separation was observed as imbalance increased. A fully paired
superfluid core formed in the trap centre, surrounded by a shell of imbalanced gas in
the normal phase.
Many properties of the zero temperature phase diagram can be deduced analytically.
On the BCS side, the two Fermi surfaces become mismatched at finite imbalance and
Cooper pairs with zero total momentum are difficult to form. When the gap between
the two Fermi surfaces becomes too large, superfluidity is broken and a first order
transition into a normal phase occurs. As the superfluid must be fully paired, the
excess atoms of the majority species will be spatially separated from the superfluid.
The result is a combination of a paired superfluid phase and a normal phase which is
a uniform non-interacting mixture of the two species, in agreement with experimental
findings. On the BEC side, the energetically favourable phase is a uniform mixture of a
superfluid gas of bosonic dimers and of a normal gas of spin polarised fermions. What
happens in the unitarity limit is not yet completely clear.
Apart from the conventional normal and superfluid phases a spin imbalance opens
the possibility for the formation of new exotic states. The most prominent one is the
FFLO phase (after Fulde, Ferrell, Larkin and Ovchinnikov). In the classical BCS picture
Cooper pairs with zero centre of mass momentum are formed by two particles at the
19
 0
 0.05
 0.1
 0.15
 0.2
 0  0.2  0.4  0.6  0.8  1
T c
/ε F
imbalance ∆n/n
superfluid
normal
Figure 1.4: Sketch of the finite temperature phase diagram of the spin imbalanced unitary
Fermi gas, as in [28]. The phase transition is first order at zero temperature
(blue square) and second order at zero imbalance (red circle), which implies the
existence of a tricritical point (green triangle). Above this point, the critical
line is second order (dashed line), and first order below (solid line).
Fermi surface. In an imbalanced system the two Fermi surfaces become mismatched.
The FFLO phase is characterised by the formation of pairs of two particles at their
respective Fermi surfaces, which consequently have non-zero centre of mass momentum.
This momentum can be interpreted as a spatial variation of the order parameter. It can
be shown that the FFLO phase occurs only in a very small domain of imbalance in the
BCS limit and is difficult to observe at unitarity. A zero-momentum pair can also be
formed if a particle on the surface of the minority Fermi sphere pairs with a particle of
opposite momentum inside the majority Fermi sphere. This phase, which is unstable,
is referred to as the Sarma phase. Exotic states due to deformed Fermi surfaces have
also been suggested.
At finite temperature thermal fluctuations reduce the gap and therefore we expect
the critical imbalance to decrease with growing temperature, resulting in a critical line
which ends in the critical point for the balanced Fermi gas. We know that at zero
temperature the phase transition is first order, whereas for the balanced gas it is second
order. This implies the existence of a tricritical point somewhere on the critical line.
A sketch of the finite temperature phase diagram for the unitary Fermi gas is given in
Fig. 1.4. The exact position of the tricritical point, as well as the shape of the critical
curve are subjects of current research. In this work we will explore the imbalanced
unitary Fermi gas close to the balanced limit and focus on the curvature of the critical
curve as the imbalance tends to zero.
20
From the numerical point of view, studying imbalance is difficult because of the
sign problem, as we shall discuss in section 4.3. The work presented here was the first
numerical study of the imbalanced unitary Fermi gas at finite temperature.
1.5 Dynamical properties
Many interesting experiments with interacting Fermi gases are performed on non-
equilibrium systems. Away from equilibrium one can access dynamical properties, for
instance the collective excitations of a gas. In particular a lot of effort has been put
into the study of spin transport phenomena in Fermi gases [33, 34, 35, 36, 37, 38]. Spin
transport plays a role in astrophysics (consider for instance neutrino transport after a
supernova explosion) and even in the physics of the early universe (hydrodynamic trans-
port of quark-gluon plasma). There is also a clear parallel to electron transport which
is fundamental for many modern technologies, like superconductors, semiconductors or
transistors, to name just a few. So far electronics has focused on charge transport, but
the electron spin might be used in a similar way as a carrier of information [39]. Either
extending conventional charge-based electronic appliances by the spin degree of free-
dom, or using the spin alone can be the foundation for a new generation of “spintronic”
devices. Advantages are for instance nonvolatility, increased data processing speed and
decreased power consumption. Understanding the spin relaxation, diffusion and other
transport properties is of fundamental importance this field.
An important advantage of cold gases in studies of spin transport phenomena is the
simplification due to the absence of relaxation mechanisms for spin currents apart from
direct collisions between atoms of different spin, unlike e.g. in a solid where collisions
with the ionic lattice can be important. The interaction strength and the temperature
of the system are also easily tunable parameters as we have discussed in section 1.1. If
the Fermi gas is strongly interacting the viscosity remains very low even well above the
critical point, which makes it easy to experimentally reach the hydrodynamic regime.
Experimental observation of different phenomena is likewise simple, as many interesting
effects happen on time scales of milliseconds and can be observed in real time. Finally,
as we shall see, very large spin polarisations can be easily created leading to large spin
currents.
In this work we will study the spin transport properties of an interacting Fermi gas
via modelling the collision of two clouds with opposite spin polarisation in an elongated
harmonic trap. This study is motivated by the recent experimental results [37, 40].
21
Chapter 2
The Fermi-Hubbard model at finite
temperature
The Fermi-Hubbard model is the simplest lattice model for two-particle scattering. It
describes non-relativistic fermions of two species, which we will label by s ∈ {↑, ↓} and
call “spin up” and “spin down”. The Hamiltonian in the grand canonical ensemble is
given by
H = H0 +H1 =
∑
k,s
(²k − µs)c†kscks + U
∑
x
c†x↑cx↑c†x↓cx↓, (2.1)
where c†ks (cks) is the time-dependent fermionic creation (annihilation) operator and
²k = 1m
∑3
j=1(1 − cos kj) the discrete dispersion relation which will be discussed in
detail in section 2.1. The first part of the Hamiltonian is the kinetic term H0 and the
second part the interaction term H1, which represents a contact (zero-range) interaction
between a spin up and a spin down particle. In the following we will assume that the
fermions have equal particle mass m and will work in units where m = 1/2.
The attractive contact interaction is characterised by the coupling constant U < 0.
This coupling can be tuned so that the scattering length takes infinite value, by solving
the two-body problem in the same way as it was done in [26]. The corresponding value
is U = −7.914 in units where m = 1/2, as will be shown in section 2.2.
We work on a 3d simple cubic spatial lattice with L3 sites and periodic boundary
conditions. The time direction remains continuous. We set the physical scale via
ν = nb3, where ν = 〈∑s c†xscxs〉 is the dimensionless filling factor, n the particle density
and b the lattice spacing. To return to the continuum, we change the chemical potentials
µs such that the gas becomes more and more dilute in lattice units, ν → 0, which then
tends to reproduce the continuum limit b → 0. Since the leading order lattice corrections
are proportional to the lattice spacing b which in turn scales like ν1/3, the continuum
extrapolation of an observable is linear with ν1/3, provided that the system is dilute
enough so that higher order corrections are negligible [26, 41]. The effect of lattice
corrections will be discussed in more depth in section 2.3. In the following we will set
the lattice spacing to unity.
22
2.1 Dispersion relation
The continuum dispersion relation for a free particle is k2/2m, corresponding to the
Fourier transform of the kinetic energy operator −∇2/2m. The according discrete
dispersion relation ²k must be a non-negative, even function in each component kj, and
due to symmetry also invariant under the permutation of the kj. In the zero-momentum
limit the discrete dispersion relation must approach the continuum one,
²k −−→k→0
k2
2m +O(k
4). (2.2)
In this work we will use the simplest lattice dispersion relation,
²k = 1m
3∑
j=1
(1− cos kj), (2.3)
for a better comparison with reference [26] where the same relation was used. This
relation can be derived directly from the following discretisation (to leading order in
lattice spacing) of the Laplace operator:
∇2cs(x) =
3∑
j=1
(cs(x+ ˆ)− 2cs(x) + cs(x− ˆ)) , (2.4)
where ˆ denote unit vectors in the three spatial directions. We can change to the
momentum space representation by expressing the fermionic fields through their discrete
Fourier transforms,
cs(x) = 1√L3
∑
p
cpseipx, (2.5)
c†s(x) =
1√
L3
∑
p
c†pse−ipx, (2.6)
where the summation goes over the first Brillouin zone, −pi < pj < pi. As a consequence
of the lattice discretisation we obtained a natural momentum cut-off. The kinetic part
23
of the Hamiltonian H0 =
∑
x c†s(x) (−∇2/2m) cs(x) then becomes
H0 = 1L3
∑
x
(
− 12m
)∑
k,p
c†kse−ikx
3∑
j=1
cps
(eip(x+ˆ) − 2eipx + eip(x−ˆ)) , (2.7)
= 1L3
∑
x
(
− 12m
)∑
k,p
c†kscpsei(p−k)x
3∑
j=1
(eipj − 2 + e−ipj) , (2.8)
= 1L3
(
− 12m
)∑
k,p
c†kscpsL3δpk
3∑
j=1
(2 cos pj − 2) , (2.9)
= 1m
∑
k
c†kscks
3∑
j=1
(1− cos kj) , (2.10)
yielding the dispersion relation (2.3). It is easy to see that as we take k → 0 this
dispersion relation tends to the correct continuum form. At large momenta however it
deviates significantly from the continuum relation. It is possible to speed the approach
to the continuum limit by choosing a more complex dispersion relation, for instance
by including higher orders in lattice spacing [42]. In section 2.3 we will discuss these
discretisation effects in more detail.
2.2 Tuning the coupling constant
The scattering amplitude for a scattering process of two non-relativistic particles at
zero temperature and chemical potential, interacting through a four-point contact in-
teraction, can be written as a series of Feynman diagrams [26, 43]:
iA(ω,p) = + . . .= + +
The Feynman rules for the diagrams are known from quantum field theory, see e.g.
[44]. Note that they depend on the metric and the definition of the propagator. In the
following we will use the (1,−1,−1,−1) metric. The Feynman rules are:
Internal lines: propagator G0(ω,p)
Vertices: factor −iU , momentum conservation
Loops: integral over undetermined loop momentum
The total incoming (and outgoing) frequency and momentum are ω and p respectively.
Since the particles are non-relativistic and at zero temperature, there is no travelling
24
backwards in time (particle number is conserved) and hence each loop carries the same
total frequency and momentum (ω,p) and has no minus sign associated with it. The
diagrammatic series can be rewritten in the following way
iA(ω,p) = +
which corresponds to
iA(ω,p) = −iU + (−iU)Π(ω,p)iA(ω,p) ⇔ A−1(ω,p) = −U−1 − iΠ(ω,p), (2.11)
where Π(ω,p) is the one-loop integral, which will be calculated below. In the low
momentum and low frequency limit the scattering amplitude is proportional to the
scattering length, as we have seen in section 1.2. Hence at unitarity A−1(0,0) = 0 and
Eq. (2.11) becomes
0 = −U−1 − iΠ(0,0) ⇒ U−1 = −iΠ(0,0). (2.12)
To obtain the coupling U we have to compute Π(ω,p), which is an integral over the
product of two propagators. To do so, we first need to calculate the free zero tem-
perature propagator G0(x − x′) = 〈0|Tc(x)c†(x′)|0〉, which we obtain from the free
non-relativistic Dirac Lagrangian,
(
i∂t + ∇
2
2m
)
G0(x− y) = iδ(4)(x− y) (2.13)
⇒ G0(ω,p) = iω − p2/2m+ iε, (2.14)
where the infinitesimal term iε specifies the contour. Accordingly, on the lattice the
propagator will have the form
G0(ω,p) = iω − ²p + iε , (2.15)
compare e.g. with [45, §9]. This explicit form depends on the metric. To calculate
the loop integral note that the loop carries an overall four-momentum (ω,p), which is
distributed over the two internal lines of the loop. Without loss of generality we can
denote these two internal momenta (ω2 + ξ, p2 + k) and (ω2 − ξ, p2 − k), so that the loop
25
integration goes over ξ and k,
Π(ω,p) =
∫ dξd3k
(2pi)4 G0(
ω
2 + ξ,
p
2 + k)G0(
ω
2 − ξ,
p
2 − k) (2.16)
=
∫ dξd3k
(2pi)4
i(
ω
2 + ξ − ²p2+k + iε
) i(
ω
2 − ξ − ²p2−k + iε
) . (2.17)
We perform the dξ integration with the residue theorem. The integrand has two simple
poles at ξ± = ±
(
ω
2 − ²p2∓k + iε
)
. We choose the counterclockwise integration in the
upper half-plane, so that the contour encloses the ξ+ pole and get
Π(ω,p) =
∫ d3k
(2pi)3
i
ω − ²p
2+k − ²p2−k
, (2.18)
dropping the iε insertion. Finally, inserting this result into Eq. (2.12) we get the
coupling constant at unitarity,
U−1 = −iΠ(0,0) = −
∫ d3k
(2pi)3
1
2²k . (2.19)
For a discrete dispersion relation the integral goes over the first Brillouin zone, −pi <
kj < pi. For the relation ²k = (1/m)
∑
j(1 − cos kj) it can be solved numerically1,
resulting in U ≈ −7.914 in units where m = 1/2.
2.3 Lattice corrections
Plugging in the expressions (2.18) and (2.19) into (2.11) we obtain an expression for
the scattering amplitude at unitarity for arbitrary frequency and momentum,
A−1(ω,p) = −U−1 − iΠ(ω,p) =
∫ d3k
(2pi)3
(
1
2²k +
1
ω − ²p
2+k − ²p2−k
)
. (2.20)
The goal is to calculate how the lattice scattering amplitude A−1(ω,p) differs from the
continuum scattering amplitude A−1cont.(ω,p) to leading order in ω and p, as this gives
an estimate of the lattice corrections. To do this we expand A−1(ω,p) via a Taylor
series around the origin, keeping only leading order terms. Introducing the following
abbreviation for the integrand,
g(ω,p) = 12²k +
1
ω − ²p
2+k − ²p2−k
, (2.21)
1One integration can be performed exactly, the remaining two numerically. The value quoted in
[26] is U ≈ −7.915.
26
we can write the first terms in the Taylor expansion as
g(ω,p) = g(0,0) + ∂g∂ω
∣∣∣∣
(0,0)
ω +
∑
i
∂g
∂pi
∣∣∣∣
(0,0)
pi +
∑
i
∂2g
∂ω∂pi
∣∣∣∣
(0,0)
ωpi
+12
∑
i,j
∂2g
∂pi∂pj
∣∣∣∣
(0,0)
pipj. (2.22)
We kept the second order momentum terms because the linear terms vanish, as will be
shown below. To calculate the derivatives of the dispersion relation we introduce the
substitution z± = p/2 ± k such that ∂/∂pi = 12∂/∂z±i. After setting p = 0 this new
variable becomes z±|p=0 = ±k and hence ²z±|p=0 = ²k due to symmetry. We will now
proceed to calculate all the terms in the expansion. The first two terms equal
g(0,0) = 12²k −
1
2²k = 0, (2.23)
∂g
∂ω
∣∣∣∣
(0,0)
ω = −1(
ω − ²p
2+k − ²p2−k
)2
∣∣∣∣∣∣∣
(0,0)
ω = −ω(2²k)2 . (2.24)
Now let us calculate the first derivative with respect to the momentum,
∂g
∂pi =
1
2
(∂²z+/∂z+i + ∂²z−/∂z−i
)
(
ω − ²p
2+k − ²p2−k
)2 . (2.25)
Since ²z± is an even function, its derivative is an odd function. Hence
∂²z+
∂z+i
+ ∂²z−∂z−i
∣∣∣∣
z±=±k
= 0. (2.26)
Therefore the contributions (∂g/∂pi)|(0,0) pi vanish and for the same reason the contri-
butions (∂2g/∂ω∂pi)|(0,0) ωpi also vanish. We have one more derivative to calculate
∂2g
∂pi∂pj =
1
4
(
ω − ²p
2+k − ²p2−k
)4
[(
ω − ²p
2+k − ²p2−k
)2( ∂2²z+
∂z+i∂z+j
+ ∂
2²z−
∂z−i∂z−j
)
+
+
(∂²z+
∂z+i
+ ∂²z−∂z−i
)
2
(
ω − ²p
2+k − ²p2−k
)(∂²z+
∂z+j
+ ∂²z−∂z−j
)]
. (2.27)
The second term in the square brackets must vanish due to (2.26). Since the second
derivative of the even function ²z± is even again, the two contributions for z+ and z−
27
are identical. Using this we obtain
∂2g
∂pi∂pj
∣∣∣∣
(0,0)
= 2(∂
2²k/∂ki∂kj)
4(2²k)2 , (2.28)
and putting everything together,
A−1(ω,p) = −
∫ d3k
(2pi)3
(
ω
(2²k)2 −
∑
i,j
(∂2²k/∂ki∂kj)pipj
4(2²k)2
)
. (2.29)
This formula holds generally for every symmetric dispersion relation and hence also for
the continuum relation k2/2m. Using this we can calculate the lattice corrections,
A−1(ω,p)−A−1cont.(ω,p) =
ω
4A−
p2
16B, (2.30)
where the values of the two parameters A and B are given below. Note that mixed
terms proportional to pipj for i 6= j make no contribution to the lattice corrections. The
reason is the symmetry of the dispersion relation: the function ∂2²k/∂ki∂kj is odd in
both ki and kj and hence its integral over the Brillouin zone vanishes. And the second
derivative ∂2/∂pi∂pj of p2/2m is only non-vanishing for i = j. Explicitly we have the
following expressions for the correction parameters:
A =
∫ d3k
(2pi)3
1
(k2/2m)2 −
∫
BZ
d3k
(2pi)3
1
²2k
, (2.31)
B =
∫ d3k
(2pi)3
1/m
(k2/2m)2 −
∫
BZ
d3k
(2pi)3
∂2²k/∂ki∂ki
²2k
, (2.32)
for arbitrary i = 1, 2, 3, c.f. [26]. In principle the leading order lattice corrections can
be completely suppressed by choosing a dispersion relation for which A = B = 0.
Such dispersion relations do indeed exist [46], but will be not explored in this work.
Comparing to the explicit expression for the effective range re in a lattice model [47]
we see that A = rem2/(2pi). Thus A = 0 is equivalent to a dispersion relation with zero
effective range. For the dispersion relation (2.3) the effective range equals
re = 8pi
∫
R\BZ
d3k
(2pi)3
1
k4 + 8pi
∫
BZ
d3k
(2pi)3
( 1
k4 −
1
(2m²k)2
)
= −0.3056, (2.33)
in units of lattice spacing [47, 48]. For this dispersion relation the effective range is
negative.
28
2.4 The finite temperature formalism
The finite temperature behaviour of a system in equilibrium is described by its partition
function Z = Tre−βH , where β is the inverse temperature. Physical observables are
thermal expectation values of operators, which are obtained through
〈Xˆ〉 = TrXˆe
−βH
Z . (2.34)
Note that the expression e−βH bears a similarity to the expression eiHt, which often
occurs in quantum field theories in relation to the time evolution operator. From
quantum field theory on the other hand we know the following identity (see e.g. [44])
for the time-evolution operator U(t), t > 0, in the interaction picture:
U(t) ≡ eiH0te−iHt = T exp
(
−i
∫ t
0
H1(t′)dt′
)
, (2.35)
where H1(t) is the interaction part of the Hamiltonian in the interaction picture, defined
by
H1(t) = eiH0tH1e−iH0t, (2.36)
and T the time-ordering operator, which places operators evaluated at later times to
the left. Equation (2.35) can be rewritten as
e−iHt = e−iH0tT exp
(
−i
∫ t
0
H1(t′)dt′
)
. (2.37)
To relate the two expressions e−iHt and e−βH we make an analytic continuation to
imaginary time, by setting all time variables t′ → −iτ ′ for a real τ ′ > 0, and identifying
τ = it with β. Doing so we can use (2.37) to write
Z = Tre−βH = Tr
[
e−βH0Tτ exp
(
−
∫ β
0
dτH1(τ)
)]
, (2.38)
where Tτ is now the imaginary time ordering operator that moves operators at larger
values of τ to the left. The exponential in equation (2.38) can be expanded in powers
of H1, as it is usually done in perturbation theory. Of course, we cannot neglect higher
order terms for a strongly coupled system, but we can use the expansion for a Monte
Carlo calculation, by randomly sampling terms in the series. Inserting the explicit form
of H1 into the expansion,
Z =
∞∑
p=0
(−U)p
∑
x1,...xp
∫
0<τ1<...β
p∏
j=1
dτjTr
[
e−βH0
p∏
j=1
c†↑(xj, τj)c↑(xj, τj)c†↓(xj, τj)c↓(xj, τj)
]
,
(2.39)
29
Z = 1 + + +− − ± . . .
Figure 2.1: Diagrammatic expansion of the partition function
generates a series of Feynman diagrams, as shown in Fig. 2.1. The diagrams consist
of four-point vertices representing a factor of (−U) and lines representing a free finite
temperature single-particle propagator [45],
Gs(0)(xi − xj, τi − τj) ≡ −〈Tτc†xis(τi)cxjs(τj)〉 (2.40)
= −Tr[Tτe−βH0c†xis(τi)cxjs(τj)], (2.41)
resulting from contractions between the fermionic operators. Lines in the upper half-
plane belong to spin up and lines in the lower half-plane to spin down propagators
respectively. Additionally, each closed fermionic loop carries a minus sign due to the
fermionic anticommutation relations. Therefore diagrams with an even and those with
an odd number of loops have different signs. This makes a direct sampling of the
individual diagrams difficult, since a series of terms with alternating signs fluctuates
and the error of the numerical calculation is bound to be large. In a certain case
however, this sign problem can be avoided, as we shall prove in the next section. The
explicit form of the propagator (2.40) is given by [49]
Gs(0)(k, τ ≡ τj − τi) =
{
e−(²k−µs)τ (1− nks) for τ > 0
−e−(²k−µs)τnks for τ ≤ 0
, (2.42)
where nks = (1 + eβ(²k−µs))−1 is the occupation of the state (k, s) for free fermions at
inverse temperature β. In the case of equal chemical potentials for spin up and spin
down, nks does not depend on s. The spatial representation of the Green’s function
can be obtained by an appropriate Fourier transform.
2.5 Matrix notation
To avoid the sign problem mentioned in the previous section, we will group the diagrams
of the expansion of Z in a clever way [50]. Consider the two parts of a Feynman diagram
corresponding to spin up and spin down separately. As the spin up and the spin down
contributions are independent each diagram can be decomposed as a product of its spin
up and spin down components. The set of all diagrams of a certain order is the sum of
all such possible products.
Now take the entity of all diagrams of arbitrary order p and consider only one
30
of the two components with spin s. Every such diagram component has p vertices
at positions x1, . . . , xp, each with an outgoing and an incoming line, which have to be
joined together somehow. We can write the set of all possible combinations of how these
lines can be connected as a p × p matrix As, where the matrix element Asij denotes a
free single-particle propagator (2.40) going from xi to xj,
Asij = Gs(0)(xi − xj, τi − τj), i, j = 1, . . . , p, (2.43)
or in other words, a connection from the outgoing line of vertex i to the incoming line of
vertex j. In this notation, each single diagram corresponds to the product of p elements
of the matrix As, with the condition that we must choose one and only one element
from each row and column of the matrix. This condition incorporates the fact that
each vertex has one and only one incoming and outgoing spin s line respectively. Such
a choice of matrix elements can be written as {Asiσ(i)|i = 1, . . . , p}, where σ ∈ Σp is
some permutation of the integers 1, . . . , p (an element of the permutation group Σp).
This permutation uniquely determines the form of the diagram.
Additionally, each loop brings a minus sign. We need a prescription how to deter-
mine the overall sign of a diagram from the permutation σ. To achieve this consider
first a random permutation and then swap two of its elements. In terms of diagrams
this means the following: If the two corresponding vertices belonged to the same spin
s loop, this loop gets broken up in two separate loops; if they belonged to two different
loops, these loops get joined together. In any case the number of loops changes by one,
which means a change in the overall sign of the diagram. Mathematically on the other
hand, such a swap means a change of the signature of the permutation. So depending
on the diagram order p the overall sign of any diagram σ of this order is either sgn(σ)
or −sgn(σ). To find out which of the two is the case, consider the identity permutation,
which corresponds to a diagram where each vertex is connected only to itself. This
means that there are p (spin s) loops. If p is even, such a diagram has sign +1, and if
p is odd, it has sign −1. In total we can write each diagram as (−1)psgn(σ)∏i Asiσ(i).
Now it is easy to see that summing all diagrams of one order we simply obtain the
determinant of the matrix As:
∑
σ∈Σp
(−1)psgn(σ)
p∏
i=1
Asiσ(i) = (−1)p detAs. (2.44)
Putting the two spin components together we get an expression for the contribution of
all diagrams of a given order p,
(−1)2p detA↑ detA↓ = detA↑ detA↓, (2.45)
where A↑ comes from the spin up part and A↓ from the spin down part of the diagrams.
31
=− +−
·
A↑11 A↑22 A↑12 A↑21
−
A↓12 A↓21A↓11 A↓22
−
Figure 2.2: Decomposition of the sum of all p = 2 diagrams into a product of two matrix
determinants.
Hence we have represented the contribution from all diagrams of a certain order p as a
product of two matrix determinants, one for each species of spin. An example of such
a representation for p = 2 is given in Fig. 2.2.
Now it is easy to see how this formalism can help with the sign problem. If the
number of spin up and spin down particles is equal (µ↑ = µ↓) then so are the Green’s
functions (2.42) for the two components and consequently also all corresponding matrix
elements. In this case we have detA↑ detA↓ = | detA|2, which is always positive. Hence
we have represented the partition function as a series of only positive terms.
Another advantage of this formalism is that we do not need to consider individual
diagrams, but obtain the contribution from all diagrams of a given order p and at a
given set of spacetime positions for the vertices (vertex configuration) at once. If we
label a vertex configuration as Sp = {(xj, τj)|j = 1, . . . , p}, the expansion of Z takes
the form
Z =
∞∑
p=0
∑
Sp
(−U)p detA↑(Sp) detA↓(Sp), (2.46)
where summation goes over all diagram orders p and vertex configurations Sp,
∑
Sp
≡
∑
x1,...,xp
∫
0<τ1<...<τp<β
p∏
j=1
dτj. (2.47)
This summation, together with the summation over p will be performed stochastically
via a Monte Carlo algorithm. To obtain physically relevant quantities we need to
compute expectation values over Z rather than Z itself. These can however be expanded
in a similar way, which will be illustrated on an example in section 4.1.
32
Chapter 3
Order parameter and finite-size scaling
To calculate the critical temperature we need to introduce an order parameter which
distinguishes between the normal and the superfluid phase. The latter is characterised
through the presence of anomalous correlations, in our case between pairs of fermions
with opposite spin, as these are the relevant degrees of freedom of the system. We
therefore first introduce the pair creation and annihilation operators
P †(x, τ) = c†x↑(τ)c†x↓(τ), (3.1)
P (x, τ) = cx↑(τ)cx↓(τ), (3.2)
and their two-point correlation function
G2(x, τ ;x′, τ ′) =
〈TτP (x, τ)P †(x′, τ ′)
〉 (3.3)
= 1ZTr[TτP (x, τ)P
†(x′, τ ′)e−βH ]. (3.4)
For an infinite system the correlation function at the critical point decays as a power
law at large separations |x− x′| → ∞,
G2(x, τ ;x′, τ ′) ∝ 1|x− x′|d−2+η , (3.5)
where d = 3 is the number of spatial dimensions and η ≈ 0.038 [51, 52] is the anomalous
dimension for the U(1) universality class. In our analysis we are limited to finite systems
characterised by the lattice size L. Finite-size effects distort the location of the critical
point by shifting the critical temperature to another “pseudocritical” value [53]. Hence
we need a prescription how to obtain Tc in the thermodynamic limit from the lattice
data.
The finite-size scaling hypothesis [53, 54] states that near the critical point the
behaviour of a finite system is determined by the dimensionless ratio L/ξ∞ between
the finite-size scale L and the correlation length of the infinite system ξ∞. From the
renormalisation group framework we know that as we approach the critical point the
correlation length diverges as ξ∞ ∝ |t|−νξ for a second order phase transition, where
t = (T − Tc)/Tc and νξ is a universal critical exponent. Hence the scale ratio L/ξ∞ is
33
 0
 0.02
 0.04
 0.06
 0.08
 0.1
 0.12
 6.5  7  7.5  8
R(
L,β
)
β
L=8
L=10
L=12
L=14
L=16
Figure 3.1: Example plot of R(L, β) as a function of the inverse lattice temperature β near
the critical point for different values of lattice size L. This data was taken at
the lattice chemical potential µ = 0.4. The lines are linear fits of the data. Due
to corrections from irrelevant operators the fits do not cross in the same point.
proportional to L|t|νξ near the critical point. The finite-size dependence will thus be
represented by a universal scaling function with the argument L|t|νξ , or equivalently by
another scaling function with the argument x = tL1/νξ .
For instance let us look at the average of G2 over space and imaginary time,
K(L, T ) = (βL3)−2
∑
x,x′
∫ β
0
dτ
∫ β
0
dτ ′G2(x, τ ;x′, τ ′), (3.6)
which is the analogue of the susceptibility [52]. Because of the power law decay (3.5) it
is clear from dimensional arguments that K(L, T ) must scale like L−(1+η). We therefore
define the rescaled averaged correlation function
R(L, T ) = L1+η(βL3)−2
∑
x,x′
∫ β
0
dτ
∫ β
0
dτ ′G2(x, τ ;x′, τ ′). (3.7)
Leaving aside corrections due to irrelevant operators for the time being, R(L, T ) near
the critical point must be a purely universal scaling function of x = tL1/νξ , analytic at
x = 0. Hence if no corrections due to irrelevant operators were present, R(L, T ) would
be independent of the lattice size L at the critical point t = 0. In this case, all R(L, T )
curves for different values of L would cross in a single point.
By looking at the example plot in Fig. 3.1 we can see that this is not the case. This
34
-3 -2 -1 1 2 3 c
-0.2
-0.1
0.1
0.2
g -gŽ
g
Figure 3.2: The relative difference (g(Li, Lj) − g˜(Li, Lj))/g(Li, Lj) as a function of c (in
lattice units), for several values of Li and Lj ranging between 10 and 16.
 0.1
 0.11
 0.12
 0.13
 0.14
 0.15
 0.16
 0  0.0005  0.001  0.0015  0.002  0.0025  0.003  0.0035  0.004
T c
g~(Li,Lj)
using L=12,14,16
using L=10,12,14,16
using all L=8,10,12,14,16
Figure 3.3: The crossing temperatures (in lattice units) for different pairs of lattice sizes
for the system from Fig. 3.1 as a function of g˜(Li, Lj). A linear dependence
is not visible. The lines represent fits through different subsets of points. For
comparison, the red box is the continuum value (with error) obtained with the
improved method, which gives consistent results independent of the lattice size
range.
35
effect is due to corrections coming from the irrelevant operators. Apart from the op-
erators present in the original Hamiltonian any renormalisation group transformation
introduces a potentially infinite hierarchy of additional operators with corresponding
coupling constants. These operators are also associated with universal critical expo-
nents under the renormalisation group flow. If the exponents are negative then the
contribution of the operator vanishes at the critical point. In this case the operator
is called irrelevant [53, 54]. Nevertheless, for finite-size systems the contribution from
these operators can become significant, for instance if the critical exponent is close to
zero. Taking the leading order correction into account, the function R(L, T ) near the
critical point can be written as a product of the universal scaling function f(tL1/νξ) and
a correction term due to irrelevant operators,
R(L, T ) = f(tL1/νξ)(1 + cL−ω + . . .), (3.8)
where c is a non-universal constant and ω another universal critical exponent. The val-
ues of the critical exponents can be determined to high precision with various methods,
see e.g. [51, 52], and for our model equal νξ ≈ 0.67 and ω ≈ 0.8. Near the critical point
we can expand Eq. (3.8) in powers of the rescaled temperature t and keeping only terms
linear in t and the leading order correction from irrelevant operators we obtain
R(L, T ) = (f0 + f1(T − Tc)L1/νξ)(1 + cL−ω). (3.9)
Previous work [26, 23] used a two-step procedure for determining Tc from this equation.
This procedure is based on the observation that although the crossings of the R(L, T )
curves do not occur exactly at Tc they are nevertheless close to the critical point for
sufficiently large L. By equating R(Li, Tij) = R(Lj, Tij) and using Eq. (3.9), a relation
between Tc and the crossing temperatures Tij can be derived,
Tij − Tc = κg(Li, Lj), (3.10)
where
g(Li, Lj) = (Lj/Li)
ω − 1
L
1
νξ
+ω
j
(
1− (Li/Lj)
1
νξ
)
+ cL
1
νξ
j
(
1− (Li/Lj)
1
νξ
−ω) , (3.11)
and κ = cf0/f1 is a non-universal constant. Note that the second term in the denom-
inator of (3.11) also contains the non-universal constant c. In the references [26] and
[23] this term is neglected, so that the second step of the procedure simplifies to a linear
fit of the crossing temperatures Tij to the function
g˜(Li, Lj) = (Lj/Li)
ω − 1
L
1
νξ
+ω
j
(
1− (Li/Lj)
1
νξ
) . (3.12)
36
Figure 3.4: A typical fit of the rescaled correlation function R(L, T ) according to Eq. (3.9).
Here data was taken at four different lattices sizes and temperatures. For this
fit χ2/d.o.f.= 1.4. The value for c was found to be −1.4(5). All quantities are
given in lattice units.
The critical temperature is the intercept of this fit. This simplification is justified if the
constant c is sufficiently small. But if c assumes values of order of unity the systematic
error associated with this approximation can reach up to 20%, as shown in Fig. 3.2. In
this case obtaining a good linear fit of the crossing temperatures is no longer possible.
An example is shown in Fig. 3.3, where the crossings from the curves depicted in Fig. 3.1
were used. No linear dependence of the crossing temperatures on g˜ is visible, moreover
depending on which values for L are included in the fit the extrapolations yield very
different values for Tc(L →∞).
To avoid this systematic uncertainty we propose a different procedure for extracting
Tc from the numerical data. In our analysis we use Eq. (3.9) directly to fit all data
triplets (R,L, T ) to a single function. An example of such a non-linear fit is shown
in Fig. 3.4, where we used the same system as for the Figures 3.1 and 3.3. The new
procedure has several advantages. Firstly, the previously described systematic error
is no longer present, which can become relevant, since we found |c| > 1 in several
cases (see Sec. 6.2 for a detailed analysis). Secondly, all information obtained from
the simulations is used for the data analysis. In the original two-step procedure the
(R, T ) tuples for each L were fitted to a line separately, which involved two unknown
parameters for each value of L. After the crossings of these lines were determined, the
information about the values of R could no longer be used for the next stage of the
analysis. From the crossings another linear fit involving two unknown parameters had
to be made. For our example from Fig. 3.4 with 16 datapoints, the original procedure
would require four independent linear fits (8 parameters) and another linear fit into
37
which the errors of the previous fits propagate. The new method suggested here only
requires a single non-linear fit of 4 parameters: f0, f1, c and Tc, of which only Tc is of
interest.
We made several consistency checks of the new method. We established that the fit
procedure gives consistent results, independent of which combinations of lattice sizes
are included in the fit, provided that all L are sufficiently large. We also confirmed that
the fit results agree if we include second order corrections in t and L−ω. For the final
results however only leading order terms were used. Also the fits are insensitive to the
precise values of the critical exponents, so that we can safely neglect the small errors
associated with their numerical determination.
The value for Tc obtained through this procedure is in lattice units and needs to
be translated into physical units. Since the only physical length scale at unitarity
is determined by the density, the corresponding physical quantity has to be Tc/εF ,
where the Fermi energy is defined as εF = (3pi2ν)2/3. In the grand canonical ensemble
the chemical potential is fixed and the corresponding filling factor ν is measured for
different values of lattice size. For sufficiently large lattices the values ν(L) scale linearly
with 1/L and an extrapolation to 1/L → 0 will yield the thermodynamic limit for the
filling factor at a given chemical potential [26]. This value will be used to determine
εF . Finally, as we discussed, the continuum limit for the critical temperature Tc/εF is
taken by extrapolating to ν → 0.
38
Chapter 4
Diagrammatic Monte Carlo sampling
4.1 Formalism
We will now introduce the mathematical formalism that will allow us to apply Monte
Carlo importance sampling techniques to the diagrammatic expansions of the partition
function and thermal expectation values. In section 2.5 we established the matrix
notation for the expansion of the partition function. The same formalism can be applied
to thermal averages of operators,
〈Xˆ〉 = 1ZTr[Xˆe
−βH ], (4.1)
which correspond to physical observables. The diagrammatic expansion of Tr[Xˆe−βH ]
will have a similar form as the one for Z, but will contain additional vertices due to
contractions with the creation and annihilation operators of which Xˆ is composed. The
general form of a diagrammatic expansion of an observable X = 〈Xˆ〉 can be written as
X = 1Z
∑
p,Sp
D(X)(Sp), (4.2)
where summation goes over all diagram orders p and vertex configurations Sp. The
diagram corresponding to the vertex configuration Sp is denoted by D(X)(Sp). Instead of
uniformly sampling the configuration space and calculating the sum, an efficient Monte
Carlo simulation will make use of importance sampling and create the configurations
Sp according to their natural probability distribution, which is given by the expansion
of the partition function:
ρ(Sp) = 1ZD
(Z)(Sp). (4.3)
Dividing by Z = ∑p,Sp D(Z)(Sp) ensures normalisation and
D(Z)(Sp) = (−U)p detA↑(Sp) detA↓(Sp), (4.4)
39
as given by equation (2.46). We can write
X = 1Z
∑
p,Sp
D(X)(Sp) =
∑
p,Sp
D(X)(Sp)
D(Z)(Sp)
D(Z)(Sp)
Z =
∑
p,Sp
D(X)(Sp)
D(Z)(Sp) ρ(Sp). (4.5)
If we create the configurations Sp according to the probability distribution ρ(Sp), we
get an estimator for the observable X through
X =
〈D(X)(Sp)
D(Z)(Sp)
〉
ρ
≡ 〈Q(X,Z)〉ρ , (4.6)
where 〈. . .〉ρ stands for averaging over a sequence of Monte Carlo vertex configurations,
created according to the probability distribution ρ(Sp) and
Q(X,Z)(Sp) = D
(X)(Sp)
D(Z)(Sp) . (4.7)
Since the terms D(X)(Sp) in the diagrammatic expansion of a physical observable are of
a similar form as D(Z)(Sp), the terms in the expansion of the partition function, we can
expect cancellations, so that the Monte Carlo estimator Q(X,Z)(Sp) will have a simple
form. To illustrate this we consider as an example the estimator for the filling factor
ν =
∑
s
〈c†xs(τ)cxs(τ)〉 =
1
Z
∑
s
Tr [c†xs(τ)cxs(τ)e−βH
] , (4.8)
where due to translational invariance (x, τ) can be an arbitrary spacetime point [26].
The diagrammatic expansion of the trace is similar to the expansion of the partition
function, the only difference being the two additional operators c†xs(τ) and cxs(τ), which
can also form contractions with the operators coming from the expansion of the Hamil-
tonian. It is easy to see that the corresponding diagrams have the form
D(ν)(Sp) =
∑
s
(−U)p detBs(Sp,x, τ) detAs(Sp), (4.9)
where s and s are opposite spin indices. The free Green’s function matrix As is the same
as in the expansion of Z and the matrix Bs is formed from the matrix As by addition
of an extra row and an extra column, originating from contractions with the operators
c†xs(τ) and cxs(τ). The Monte Carlo estimator of the number density is consequently
given by the ratio of two determinants,
Q(ν,Z)(Sp) = D
(ν)(Sp)
D(Z)(Sp) =
∑
s
detBs(Sp,x, τ)
detAs(Sp) . (4.10)
40
Monte Carlo estimators for other thermodynamic observables can be obtained in an
analogous way.
4.2 Extending the configuration space
As we will elaborate in the next chapter, one of the core ideas of the algorithm is the
extension of the space of allowed vertex configurations Sp. Mathematically, this means
extending the domain on which the sampling probability distribution is defined. In
this section we will develop the corresponding mathematical formalism. Let us call the
initial domain S(Z) and the extended domain S(ZW ) = S(Z) ∪ S(G). Then we define
the extended probability distribution to be
D(ZW )(Sp) =
{
D(Z)(Sp) for Sp ∈ S(Z)
D(G)(Sp) for Sp ∈ S(G)
(4.11)
in terms of the diagrammatic expansion. Apart from normalisation the probability
distribution stays the same on the initial domain. Taking normalisation into account:
ρW (Sp) ≡ 1ZW D
(ZW )(Sp) =
{
ρ(Sp)Z/ZW for Sp ∈ S(Z)
ρG(Sp) for Sp ∈ S(G)
, (4.12)
where Z = ∑S(Z) D(Z)(Sp) and ZW =
∑
S(ZW ) D(ZW )(Sp). In other words, the new
probability distribution on S(Z) is obtained from the original one by reweighting with
the ratio of the measures of the old and the new domains.
Thermal averages must stay invariant under this extension, to ensure that the phys-
ical predictions are still valid. If the observable X is defined on S(Z) only, we extend it
to S(ZW ) by setting D(X)(Sp) = 0 ∀Sp ∈ S(G). This implies
D(X)(Sp)
D(Z)(Sp) =
D(X)(Sp)
D(ZW )(Sp) (4.13)
⇒ Q(X,Z)(Sp) = Q(X,ZW )(Sp) ∀Sp ∈ S(ZW ). (4.14)
Applying (4.12) and (4.13) to (4.5) we can estimate X through
X =
∑
p,Sp
D(X)(Sp)
D(ZW )(Sp)ρW (Sp)
ZW
Z =
ZW
Z 〈Q
(X,ZW )〉ρW . (4.15)
As far as the partition function Z is concerned,
Z =
∑
p,Sp
D(Z)(Sp) = ZW
∑
p,Sp
D(Z)(Sp)
D(ZW )(Sp)
D(ZW )(Sp)
ZW = ZW 〈Q
(Z,ZW )〉ρW , (4.16)
41
and hence Z
ZW = 〈Q
(Z,ZW )〉ρW . (4.17)
But since
D(Z)(Sp) =
{
D(ZW )(Sp) for Sp ∈ S(Z)
0 for Sp ∈ S(G)
(4.18)
we obtain a very simple form for the corresponding Monte Carlo estimator,
Q(Z,ZW )(Sp) =
{
1 for Sp ∈ S(Z)
0 for Sp ∈ S(G)
(4.19)
Now we have all the necessary preliminaries to estimate physical observables on the
extended domain:
X = 〈Q(X,Z)〉ρ = ZWZ 〈Q
(X,ZW )〉ρW =
〈Q(X,ZW )〉ρW
〈Q(Z,ZW )〉ρW
. (4.20)
4.3 Sampling with a different distribution
The original DDMC algorithm relies strongly on the assumption of equal densities of the
two fermionic species. This assumption allows us to write the partition function (2.46)
as a sum of positive terms only, and consequently to use it as a probability distribution
for Monte Carlo sampling. In the presence of an imbalance µ↑ 6= µ↓ a generalisation
of the algorithm is necessary, as the diagrams D(Z)(Sp) no longer need to be positive
and hence ρ(Sp) ceases to be a well-defined probability distribution. This forces us to
use a different distribution for Monte Carlo sampling. The mathematical formalism
for sampling with respect to a different probability distribution ρ′(Sp) = 1Z′D(Z
′)(Sp),
defined on the same domain as ρ(Sp), will be introduced in this section. Proceeding as
in (4.5) we can write
X = Z
′
Z
∑
p,Sp
D(X)(Sp)
D(Z′)(Sp)ρ
′(Sp) = Z
′
Z 〈Q
(X,Z′)〉ρ′ . (4.21)
To obtain an expression for the prefactor consider the diagrammatic expansion of Z in
terms of the new distribution ρ′(Sp)
Z =
∑
p,Sp
D(Z)(Sp) = Z ′
∑
p,Sp
D(Z)(Sp)
D(Z′)(Sp)ρ
′(Sp) ⇒ ZZ ′ = 〈Q
(Z,Z′)〉ρ′ . (4.22)
Consequently
X = 〈Q
(X,Z′)〉ρ′
〈Q(Z,Z′)〉ρ′
, (4.23)
42
which is analogous to (4.20). Finally, as discussed in section 4.2, we can extend the
modified probability distribution ρ′(Sp) to ρ′W (Sp) on a larger domain, in the same way
as we did for ρ(Sp). We obtain using (4.20),
〈Q(X,Z′)〉ρ′ =
〈Q(X,Z′W )〉ρ′W
〈Q(Z′,Z′W )〉ρ′W
(4.24)
and
〈Q(Z,Z′)〉ρ′ =
〈Q(Z,Z′W )〉ρ′W
〈Q(Z′,Z′W )〉ρ′W
, (4.25)
and inserting everything into equation (4.23),
X =
〈Q(X,Z′W )〉ρ′W〈Q(Z′,Z′W )〉ρ′W
〈Q(Z′,Z′W )〉ρ′W〈Q(Z,Z′W )〉ρ′W
=
〈Q(X,Z′W )〉ρ′W〈Q(Z,Z′W )〉ρ′W
. (4.26)
To avoid the sign problem an obvious choice for an alternative distribution is the “sign
quenched” distribution. The idea is based on the “phase quenched method” known
from lattice QCD [55]. We can write the function ρW (Sp) as a product of its modulus
and its sign,
ρW (Sp) = 1ZW (−U)
p| detA↑(Sp) detA↓(Sp)|sgn(Sp), (4.27)
and use the positive function
ρ′W (Sp) ≡
1
Z ′W
(−U)p| detA↑(Sp) detA↓(Sp)| (4.28)
as the new probability distribution. Using equation (4.26) we can immediately read off
the general expression for a thermal average of an operator Xˆ:
〈Xˆ〉ρ =
〈Q(X,Z′W )〉ρ′W〈Q(Z,Z′W )〉ρ′W
=
〈Q(X,ZW )(Sp)sgn(Sp)
〉
ρ′W
〈Q(Z,ZW )(Sp)sgn(Sp)〉ρ′W
. (4.29)
The Monte Carlo estimators Q remain unchanged, apart from a multiplication with
±1 depending on the relative sign of the two matrix determinants detA↑(Sp) and
detA↓(Sp).
43
Chapter 5
Implementation of the algorithm
Now that we have introduced the diagrammatic Monte Carlo formalism and showed
how Monte Carlo estimators for physical observables can be obtained, we are ready to
present the algorithm through which the sampling is realised. The configuration space
is sampled via a Monte Carlo Markov chain process: in each step one of the possible
updates from the current vertex configuration Sp to another vertex configuration S ′q is
proposed with probability W (Sp → S ′q), and accepted with probability P (Sp → S ′q) =
min(1,R), given by the detailed balance equation,
RW (Sp → S ′q)D(Z
′
W )(Sp) = W (Sp ← S ′q)D(Z
′
W )(S ′q). (5.1)
The requirements of detailed balance and ergodicity (from any initial configuration
every other vertex configuration can be reached in a finite number of steps) ensure that
the produced configurations are indeed distributed according to the correct thermal
probability distribution ρ′W (Sp).
The diagrammatic expansion of the two-point correlation function (3.4) is similar
to the expansion of the partition function Z, but contains an additional pair of two-
point vertices at (x, τ) and (x′, τ ′). It is thus of advantage to sample these two series
in the same simulation. In addition to sampling the regular four-point diagrams we
allow updates that insert the pair of two-point vertices (which we will call “worm
vertices”) into the configuration space. The vertex at (x, τ) contains two incoming lines
corresponding to spin up and spin down, and the vertex at (x′, τ ′) two outgoing lines
corresponding to spin up and spin down. In matrix notation this means that if the worm
vertices are present, the Green’s function matrices As each get an additional row and
column, coming from contractions with P † and P respectively. We have already worked
out the formalism for extending the configuration space in section 4.2. The original
configuration space of four-point vertices will be denoted S(Z) and the configuration
space containing the pair of two-point vertices S(G), such that S(ZW ) = S(Z) ∪ S(G).
The extended partition function ZW takes the form
ZW = Z
(
1 + ζ
∑
x,x′
∫ β
0
dτ
∫ β
0
dτ ′G2(x, τ ;x′, τ ′)
)
, (5.2)
44
 40
 60
 80
 100
 120
 140
 160
 0  20000  40000  60000  80000  100000  40
 60
 80
 100
 120
 140
 160
 0  20000  40000  60000  80000  100000
Figure 5.1: The first 100000 numerical measurements of the interaction energy (a measure-
ment takes place every 100 Monte Carlo steps) with the worm setup (left) and
the diagonal setup (right). Strong autocorrelations are visible in the worm setup.
 0
 0.002
 0.004
 0.006
 0.008
 0.01
 0.012
 0.014
 0  10000  20000  30000  40000  50000  60000  70000
Re
lati
ve 
Err
or
Block Size
 0
 0.002
 0.004
 0.006
 0.008
 0.01
 0.012
 0.014
 0  10000  20000  30000  40000  50000  60000  70000
Re
lati
ve 
Err
or
Block Size
Figure 5.2: The blocking error analysis of the interaction energy with the worm setup (left)
and the diagonal setup (right). The blocked error is much higher in the worm
setup and continues increasing even for large block sizes.
where ζ is an arbitrary parameter. The advantage of this setup is that the Monte Carlo
estimator for the order parameter R(L, T ) defined in equation (3.7) now becomes very
simple: it is just a constant times the ratio of configurations with and without worm
vertices,
Q(R,ZW )(Sp) =
{
0 for Sp ∈ S(Z)
L1+η(βL3)−2ζ−1 for Sp ∈ S(G)
. (5.3)
Other physical observables like the number density or the energy are still only measured
when the system is in the “physical sector” S(Z), namely when the worm vertices are
not present.
A detailed description of the original worm updates can be found in the appendix of
[26]. The main idea behind the algorithm is that at low densities the major contribution
to the partition function comes from multi-ladder diagrams – these are configurations
where the vertices are arranged into several vertex chains. Proposing updates that
favour the creation of such vertex chains will lead to higher acceptance ratios and
45
increase the efficiency of the simulation. At low densities the acceptance ratios of the
worm updates are an order of magnitude higher than the acceptance ratios of the simple
diagonal updates, in which the vertices are inserted or removed at random.
We found however, that the worm type addition and removal updates from the orig-
inal setup [26] suffer from strong autocorrelations, so that even after many successful
updates the configuration does not change significantly [5, 6]. The cause of this problem
is that all updates are performed exclusively through the worm vertices. More precisely,
a four-point vertex can only be added or removed from the configuration if its coordi-
nates are in the vicinity of one of the worm vertices. This implies that over long times
(until the position of the worm vertices changes significantly through another update)
only a small proportion of the configuration space can be changed. Moreover the algo-
rithm tends to repeatedly accept and then immediately undo a change in consecutive
Monte Carlo steps.
To illustrate the extent of the problem we compare the measurements of the inter-
action energy (which is proportional to the diagram order) in the worm setup and the
diagonal setup in Fig. 5.1. Both simulations used the same parameters and a compara-
ble number of Monte Carlo steps. Figure 5.2 shows the blocking analysis of the relative
error for the same quantity. Blocking is a widely used technique to estimate the error
of an autocorrelated measurement. The single data points are arranged consecutively
into blocks of equal size, and each block is replaced by the average of the measurements
it contains. Then the error is calculated for the resulting blocked system in the usual
way. If no autocorrelations are present, the error will be independent of the block size.
In the presence of autocorrelations N consecutive data points fluctuate less than N in-
dependent measurements. Hence the error will increase with block size, until the block
size reaches the autocorrelation length of the system.
Because of the large errors due to autocorrelations the worm setup is effectively less
efficient than the standard diagonal setup. For this reason in the present study we em-
ploy the conventional diagonal updates, together with the modified worm addition and
removal updates, as proposed in [5, 6]. This setup combines the advantages of the di-
agonal setup (weak autocorrelations) with the ones of the worm setup (high acceptance
ratios). Below is a summary of all updates used in our simulation. For the modified
updates we also give the values of the acceptance ratios R (the other acceptance ratios
can be found in [26], in the presence of an imbalance one merely needs to replace the
terms | detA|2 by | detA↑ detA↓|).
Updates only concerning the worm vertices:
• Worm creation/annihilation: insert/remove the pair P (x, τ), P †(x′, τ ′) into/
from the configuration. In our setup the distributions for the coordinates of P
and P † are independent: both are distributed uniformly over the lattice, so that
46
W (Sp → S˜p) = (βL3)−2, where S˜p stands for the configuration Sp with the addi-
tional two-point vertices. The authors of [26] describe a setup in which the vertex
P is selected randomly, and the vertex P † is then chosen in a spacetime hypercube
of given extent around P . To avoid autocorrelations that can be associated with
this scheme we employ the independent setup. The acceptance ratio is then
R =
∣∣∣∣∣
detA↑(S˜p) detA↓(S˜p)
detA↑(Sp) detA↓(Sp)
∣∣∣∣∣ (βL
3)2ζ. (5.4)
• Worm shift: shift the P †(x′, τ ′) vertex to other coordinates. This update is
equivalent to the worm shift update in [26] and involves a shift to a nearest neigh-
bour on the lattice and a time shift in some interval around the old coordinates.
Updates of the regular four-point vertices: adding/removing a four-point vertex
(changes the diagram order).
• Diagonal version: add or remove a random vertex. This is the most basic setup
for changing the diagram order, however at low densities the acceptance ratios
are very low.
• Modified worm-type updates:
– Choose a random four-point vertex from the configuration (which will act as
a worm for this step).
– Addition: add another four-point vertex on the same lattice site and in some
time interval of length ∆τ around the worm.
– Removal: remove the nearest neighbour of the worm vertex (implies that
addition can only be accepted if the new vertex is the nearest neighbour of
the worm).
The probability density for the addition update is then W (Sp → Sp+1) = 1/(p∆τ),
where 1/p comes from selecting the worm and 1/∆τ from choosing the new time
coordinate. Analogously for the removal update W (Sp ← Sp+1) = 1/(p + 1) and
the acceptance ratio becomes
R =
∣∣∣∣
detA↑(Sp+1) detA↓(Sp+1)
detA↑(Sp) detA↓(Sp)
∣∣∣∣
(−U)p∆τ
p+ 1 . (5.5)
The modified worm setup still prolongs existing vertex chains like the original worm
setup, but autocorrelations are significantly reduced since the “worm” changes with
every update. This new type of updates can only be employed in addition to the regu-
lar diagonal addition and removal updates. It works regardless if the pair of two-point
47
 0
 0.001
 0.002
 0.003
 0.004
 0.005
 0.006
 0.007
 0.008
 0  20  40  60  80  100  120  140
Re
la
tiv
e 
Er
ro
r
Block Size
Figure 5.3: The blocking error analysis of the interaction energy. Red circles correspond
to the pure diagonal setup and blue squares to a combination of diagonal and
modified worm updates with equal probabilities for each kind of update.
vertices is present or not in the configuration (the original worm addition and removal
updates can only take place when the two-point vertices are present). The acceptance
rates for this update are comparable with those for the regular worm updates. To
demonstrate the increase in efficiency we compare the diagonal and the modified worm
setup at low density, when the acceptance rates of the diagonal updates are particularly
poor. We again consider the blocking analysis of the relative error for the interaction
energy. As Fig. 5.3 clearly shows, the autocorrelation length does not increase in pres-
ence of the modified worm updates. The blocked error in this case is significantly lower
due to the increased acceptance rate.
The structure of the updates requires the calculation of a matrix determinant of a
large matrix (with rank p up to about 6000) in each Monte Carlo step. Since only a few
elements of the matrix change in the course of one update (at most one row and one
column) a recalculation of the whole determinant from scratch is not necessary. In our
implementation we make use of the fast matrix update formulae [50] which decrease the
number of required operations to order p2 instead of order p3. A detailed description
of the matrix update procedures with all necessary mathematical proofs is given in
the appendix A. In the presence of an imbalance we need to keep in memory two large
matrices instead of one and update each of these matrices separately. Another drawback
is that the relative error of the sign adds to the relative error of each observable.
Numerical errors can become very large if the expectation value of the sign in the
denominator of a generic thermal expectation value (4.29) is close to zero, as it happens
48
0.00 0.05 0.10 0.15 0.20
DΜ
¶F
0.2
0.4
0.6
0.8
1.0
<sign>
Figure 5.4: Schematic plot of the average sign near the critical point as a function of imbal-
ance. The shaded area covers the range of values the sign can take at different
values of lattice size and chemical potential. The lower boundary of this area
is the “worst-case” curve of the sign, corresponding to the largest lattice sizes
used.
for the expectation value of the phase in QCD. For the unitary Fermi gas the sign
remains very close to unity for small imbalances, as shown in Fig. 5.4, so that sign
quenching is applicable for imbalances up to approximately 0.2εF . The restricting factor
that keeps us from reaching large imbalances is not the sign, but rather the fact that
even large values of the lattice chemical potential difference ∆µ do not necessarily lead
to large differences in the filling factors of the two components and hence the physical
value ∆µ/εF still remains small. This method works best close to the balanced limit
and can provide a useful tool to examine the trend of the critical temperature and other
observables for small deviations from it.
49
Chapter 6
Observables at the critical point
The results presented in this chapter were published in [6, 7].
With the DDMC algorithm described in the previous chapters we obtained data
at 25 different values (µ↑, µ↓), of which 8 were at µ↑ = µ↓. The lattice sizes varied
between 43 for the highest filling factor and 263 for the lowest, so that the volume
range in physical units was approximately constant. As discussed in chapter 2, the
dimensionless physical observables scale linearly with ν1/3 if we are close enough to the
continuum limit, such that higher order lattice corrections can be neglected. With our
data this behaviour is seen for ν1/3 / 0.75. This condition was fulfilled for 23 out of
the 25 points and in particular for 7 out of the 8 balanced points.
The two most common ways of quantifying imbalance are either through the chemi-
cal potential difference ∆µ/εF = |µ↑−µ↓|/εF , or through the relative density difference
∆ν/ν = |ν↑ − ν↓|/(ν↑ + ν↓). For the values of imbalance considered in our study these
two quantities are proportional to each other, with ∆ν/ν = 0.122(2)∆µ/εF , as illus-
trated in Fig. 6.1. The relative density difference shows no dependence on lattice size
(the L-dependencies of ν and ∆ν cancel each other out), but considerable dependence
on the temperature. Also since ∆ν is a very small quantity, numerical fluctuations can
become significant. Since the chemical potential difference is less prone to numerical
errors, we will use it from now on to quantify imbalance.
Every observable X is a function of filling factor ν and imbalance h = ∆µ/εF .
We are ultimately interested in the continuum limit ν = 0 and want to perform the
corresponding extrapolation. To achieve this all numerical data is fitted to a three
dimensional surface, where the following assumptions are made for the form of the
fitted function:
• At fixed imbalance, X is a linear function of ν1/3 with slope α(X)(h): X(ν, h) =
X(h) + α(X)(h)ν1/3.
• X(h) and α(X)(h) viewed as functions of the imbalance h are analytic and can
thus be Taylor expanded.
• Due to symmetry in h all odd powers in the Taylor expansions of X(h) and α(X)(h)
have to vanish.
50
 0
 0.005
 0.01
 0.015
 0.02
 0.025
 0.03
 0.035
 0  0.05  0.1  0.15  0.2  0.25
∆ν
/ν
∆µ/εF
Figure 6.1: Relation between the chemical potential difference and the relative density dif-
ference at Tc.
If for instance we expand X(h) and α(X)(h) to leading order in h the fitted function
becomes
X(ν, h) = X0 +X2h2 + (α(X)0 + α(X)2 h2)ν1/3. (6.1)
This requires a linear fit of four parameters.
6.1 Spin susceptibility and pseudogap
From the fit shown in Fig. 6.1 we can extract the spin susceptibility of the gas [56].
The dimensionless spin susceptibility χ˜ (in units of the spin susceptibility of an ideal
Fermi gas) can be obtained from the relative density difference via
∆ν
ν =
3
2 χ˜
∆µ
2εF . (6.2)
From our data we obtain the value χ˜ = 0.163(3) at the critical point. This value is lower
than the zero-temperature result χ˜(T = 0) ≈ 0.54 extracted from the experimental data
in [56].
Experimentally the spin imbalance ∆ν/ν is increased by applying an external mag-
netic field h = ∆µ/εF . At zero temperature the system does not react immediately
towards a change in h, resulting in a gap ∆. The gap can be understood as half of the
energy needed to break one fermionic pair in order to polarise the system by one excess
particle [3]. The size of the gap can be extracted from the density of states which is zero
51
in the interval [µ−∆, µ+∆]. At finite temperature the density of states can exhibit a
dip which is similar to the gap and is referred to as the pseudogap ∆∗ (see e.g. [57] for
a review). The existence of the pseudogap can also be verified by looking at the spin
susceptibility of the system. In a degenerate ideal gas which is described by the Fermi-
Dirac distribution, the susceptibility equals the density of states at the Fermi level [3].
In this simple model, since for small imbalances the Fermi levels of the two spin states
will lie within the dip, the spin susceptibility will be lower than for large imbalances,
when the Fermi levels of the two spin states are both outside the dip. As a consequence
∆ν/ν as a function of h will be a curve with two asymptotes: as the imbalance goes
to zero the asymptote has zero intercept and a small slope, while in the limit of large
imbalance the asymptote has a greater slope and the intercept is unequal zero. Our
data suggests an intercept which is either zero or very small. For a fit where both the
slope and the intercept are free parameters we obtain ∆ν/ν = 0.149(6)h − 0.0016(3).
From this we can conclude that if the pseudogap exists for the unitary Fermi gas at Tc
it must be either very broad (∆∗ & 0.2εF , the maximal imbalance considered in this
work) or the dip in the density of states is very small, such that it cannot be resolved
within the precision of our data.
6.2 Critical temperature
Before we include the imbalanced data we first present our analysis of the data at
zero imbalance, for comparison with previous results from [26]. Our results and the
continuum extrapolation are shown in Fig. 6.2. A line was fitted through the seven
points with ν1/3 < 0.75, resulting in Tc/εF = 0.173(6) − 0.16(1)ν1/3. The goodness
of fit is χ2/d.o.f. = 0.39. For comparison we also fit a quadratic through all eight
data points, resulting in a continuum value of Tc/εF = 0.188(15), which is in excellent
agreement with the linear extrapolation. This confirms that sub-leading corrections
proportional to ν2/3 can indeed be neglected for sufficiently small ν.
In Fig. 6.3 we show the results for the fit parameters c, f0 and f1Tc, according to
Eq. (3.9). These parameters are smooth functions of the filling factor. In particular
this data confirms that the non-universal constant c does indeed take values of order
unity (cf. section 3) and thus cannot be neglected.
Now we also include data with µ↑ 6= µ↓. The best fit of the critical temperature in
units of εF to the function
Tc(ν, h) = T0 + T2h2 + (α(T )0 + α(T )2 h2)ν1/3, (6.3)
yields T0 = 0.171(5), α(T )0 = −0.154(9), T2 = 0.4 ± 0.9 and α(T )2 = −0.7 ± 1.9 with
χ2/d.o.f.= 0.43. Note that the T2 value corresponding to the minimal χ2 is positive,
which is forbidden by physical arguments (Tc can only decrease with growing imbalance,
52
 0.05
 0.1
 0.15
 0.2
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9
T c
/ε F
ν1/3
Tc(ν = 0) = 0.173(6)εF
Figure 6.2: The critical temperature versus filling factor for different values of the chemical
potential. The continuum limit corresponds to ν → 0. The linear extrapolation
(solid line) of the seven data points at lowest filling factors (filled circles) yields
a continuum value of Tc/εF = 0.173(6). The dashed line corresponds to a
quadratic fit through all data points.
-6
-5
-4
-3
-2
-1
 0
 1
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9
c
ν1/3
-0.2
-0.1
 0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9
ν1/3
f(0)f’(0)
Figure 6.3: The non-universal constant c (left) and the first two coefficients in the expansion
of the universal scaling function f(x) (right) in lattice units versus filling factor,
see Eq. (3.9) with f(0) = f0 and f ′(0) = f1Tc.
53
 0.05
 0.075
 0.1
 0.125
 0.15
 0.175
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
T c
/ε F
ν1/3
Figure 6.4: Projection of the data onto the (ν1/3-Tc) plane. Red circles denote the balanced
data and blue triangles data at non-zero imbalance. The line corresponds to the
constant fit (6.7). Dashed lines denote the error margins.
as interactions are suppressed). The χ2-function is very flat along the T2 direction, so
that forcing T2 = 0 results in χ2/d.o.f.= 0.44. From the error on T2 we derive the lower
bound T2 > −0.5. The best fit values for T0 and α(T )0 are in excellent agreement with
the ones obtained from the fit of the balanced data only.
The error on the best fit value for α2 is very large and the fit is consistent with
α(T )2 = 0. Hence we also perform a fit to the function
Tc(ν, h) = T0 + T2h2 + α(T )0 ν1/3, (6.4)
where Tc(h) has again been expanded to quadratic order and the function α(T )(h) has
been replaced by a constant α(T )0 . The best fit is
Tc(ν, h) = 0.171(5) + 0.07(11)h2 − 0.155(8)ν1/3, (6.5)
with χ2/d.o.f.= 0.41. This χ2-value is even lower than for the previous fit, which means
that the data justifies dropping the α(T )2 term. The best fit result is still consistent with
T2 = 0 and leads to a much tighter lower bound T2 > −0.04. The other parameters T0
and α(T )0 agree with the results from the previous fit and the fit of the balanced data.
Since our results indicate that Tc remains almost unchanged in response to a weak
54
 0.13
 0.14
 0.15
 0.16
 0.17
 0  0.05  0.1  0.15  0.2  0.25
∆µ/εF
Tc/εF
Figure 6.5: The continuum limit of the critical temperature as a function of imbalance. The
solid line is the value obtained from the constant fit (6.6), the shaded area corre-
sponds to one standard deviation. The dashed line is the lower bound obtained
from fit (6.3) and the dot-dashed line is the tighter lower bound obtained from
fit (6.4).
Figure 6.6: Stereogram of the critical temperature versus filling factor and imbalance. The
surface corresponds to the constant fit (6.7). To view the image focus your eyes
on a point between the image and you (cross viewing).
55
T0 δT0 T2 T2 lower bound α(T )0 δα(T )0 α(T )2 δα(T )2 χ2/d.o.f.
balanced data 0.173 0.006 -0.16 0.01 0.39
fit to Eq. (6.3) 0.171 0.005 0.4 -0.5 -0.154 0.009 -0.7 1.9 0.43
fit to Eq. (6.4) 0.171 0.005 0.07 -0.04 -0.155 0.008 0.41
fit to Eq. (6.6) 0.172 0.0045 -0.156 0.008 0.41
Table 6.1: Comparison of fit parameters obtained by different fit methods.
imbalance, we also perform a fit to constant Tc(h) and α(T )(h),
Tc(ν, h) = T0 + α(T )0 ν1/3. (6.6)
This is the same function as the one used in the balanced case and corresponds to a
straight line fitted through the projection of all data points onto the (ν1/3-Tc) plane,
see Fig. 6.4. The best fit is
Tc(ν, h) = 0.1720(45)− 0.156(8)ν1/3, (6.7)
with χ2/d.o.f.= 0.41. Again the result agrees with the previous fits.
We also performed fits using the jackknife method and several robust fits. All results
were consistent with the minimal χ2 fits. Table 6.1 provides an overview of the results
obtained with the different fit methods. The values for the parameters T0 and α(T )0
were obtained with high accuracy and are all in excellent agreement with each other,
independently of the form of the fit function. Depending on the model assumptions two
lower bounds could be derived for the leading order deviation of the critical temperature
from its balanced value. Figure 6.5 shows these two bounds compared with the value
in the balanced case. A three dimensional stereoscopic plot of the data together with a
constant surface fit is presented in Fig. 6.6.
6.3 Energy per particle
The total energy is composed of the kinetic energy Ekin and the interaction energy Eint.
An explicit expression for the former can be obtained from the coordinate space picture,
Ekin = −
〈∑
x,s
c†xs∇2cxs
〉
(6.8)
= −
〈∑
x,s
c†xs
3∑
j=1
(c(x+jˆ)s + c(x−jˆ)s − 2cxs)
〉
(6.9)
=
∑
s
〈L3(6c†xscxs − 6c†xsc(x+jˆ)s)〉, for any jˆ, x. (6.10)
56
The factor L3 in the last step comes from summation over all lattice sites and the factors
of 6 are due to summation over j. Since ∑s〈c†xscxs〉 corresponds to the filling factor,
the kinetic energy per particle can be written as
Ekin/L3ν = 6
(
1−
∑
s〈c†xsc(x+jˆ)s〉
ν
)
. (6.11)
The Monte Carlo estimator for the interaction energy Eint = 〈H1〉 can be obtained with
the Hellmann-Feynman theorem [26]. Define H(λ) = H0 + λH1 and observe that
TrH1e−βH = − 1β
∂
∂λTre
−βH(λ)
∣∣∣∣
λ=1
≡ − 1β
∂Z(λ)
∂λ
∣∣∣∣
λ=1
. (6.12)
The modified function Z(λ) differs from the partition function only by a rescaling of
the coupling constant U → λU . From equation (2.39) it follows immediately that a
diagram of order p is proportional to λp and hence differentiating and then setting λ = 1
generates a factor of p. Hence the Monte Carlo estimator for the interaction energy is
given by Q(H1,Z)(Sp) = −β−1p.
The energy per particle thus depends on three observables, the kinetic correlator∑
s〈c†xsc(x+jˆ)s〉, the filling factor
∑
s〈c†xscxs〉 and the average diagram order. Naturally,
these observables are strongly correlated with each other. Calculating their expectation
values and errors separately and then using error propagation would lead to a tremen-
dous overestimation of the total error. We therefore apply a different strategy. We
group the data into N blocks of consecutive data points, with block size well above the
autocorrelation size. For each block we then calculate the average of each of the three
observables and then from these averages the value of the total energy per particle in
units of the Fermi energy. This leaves us we a set of N independent values for the energy
per particle. We then calculate the total average, which will be our best estimate for the
energy per particle and obtain the error on this estimate with the jackknife method, see
e.g. [58], as follows. We remove one block at a time from the total sample and calculate
the averages of the N reduced samples with N − 1 points each. The variance of this
sample of reduced averages yields an error estimate on the total error of the energy per
particle.
Our results using only balanced data are shown in Fig. 6.7. These results were
obtained at Tc, but the temperature dependence of the energy per particle near Tc was
found to be very weak. Also this quantity showed no dependence on the lattice size L
and hence the final results are averages over all values of L used. For the energy per
particle we obtain the continuum value E/NεF = 0.276(14). In units of the ground
state energy of the free gas, EFG = (3/5)NεF , our result is E/EFG = 0.46(2). Since
the expression for the kinetic energy (6.11) involves a difference of two quantities of
comparable size, large fluctuations can occur, especially at low filling factors. For this
57
 0.25
 0.3
 0.35
 0.4
 0.45
 0.5
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9
E/
E F
G
ν1/3
E(ν = 0) = 0.46(2)EFG
Figure 6.7: The energy per particle versus filling factor at the critical point. The linear fit
was performed for data at filling factors ν1/3 < 0.75. Filled circles indicate data
included in the fit and empty circles stand for data excluded from the fit, see
text for discussion.
reason measurements of the energy per particle at lowest filling factor could only be
performed on lattices with size L ≤ 14, which is smaller than the lattice sizes used
for the measurement of the critical temperature. We include this point in the plot in
Fig. 6.7, but exclude it from the linear fit. The goodness of fit is χ2/d.o.f. = 2.1.
Since with increasing imbalance interactions become suppressed, we expect the ab-
solute value of the interaction energy to decrease. This in turn means an increase of
the total energy, since the interaction energy is negative. As we did for the critical
temperature we fit the energy in units of EFG to the function
E(ν, h) = E0 + E2h2 + (α(E)0 + α(E)2 h2)ν1/3 (6.13)
and obtain the best fit parameters E0 = 0.440(15), α(E)0 = −0.17(3), E2 = 3.4±2.2 and
α(E)2 = −3.1± 4.5, with χ2/d.o.f.= 2.8. These results are consistent with the balanced
fit. The leading coefficient E2 = 3.4 ± 2.2 is no longer consistent with zero. We also
perform a fit to the function
E(ν, h) = E0 + E2h2 + α(E)0 ν1/3 (6.14)
and obtain the best fit result
E(ν, h) = 0.444(13) + 1.9(3)h2 − 0.18(2)ν1/3, (6.15)
58
 0.25
 0.3
 0.35
 0.4
 0.45
 0.5
 0.55
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
E/
E F
G
ν1/3
Figure 6.8: Projection of the data for the energy per particle onto the (ν1/3-E) plane. Red
circles denote the balanced data and blue triangles data at non-zero imbalance.
For comparison the fit at zero imbalance is shown. The points at non-zero
imbalance tend to lie above the balanced fit line.
Figure 6.9: Stereogram of the energy per particle E/EFG versus filling factor and imbalance.
The surface fit corresponds to the quadratic fit (6.13). To view the image focus
your eyes on a point between the image and you (cross viewing).
which agrees with the previous result. The χ2/d.o.f.= 2.7. The figures 6.8 and 6.9 show
the numerical data.
59
 0.2
 0.25
 0.3
 0.35
 0.4
 0.45
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
µ/
ε F
ν1/3
Figure 6.10: Projection of the data for the average chemical potential onto the (ν1/3-µ)
plane. Red circles denote the balanced data and blue triangles data at non-
zero imbalance. The solid line corresponds to the constant fit and the dashed
lines indicate the error margins.
6.4 Chemical potential
For the chemical potential at Tc we obtain the continuum value µ/εF = 0.429(9) with
χ2/d.o.f. = 2.8 using only balanced data. A similar analysis can be performed for the
average chemical potential µ/εF = |µ↑+µ↓|/2εF in the presence of an imbalance. Since
the average chemical potential is not expected to depend on the imbalance we fit our
data in units of εF to the constant function
µ(ν, h) = µ0 + α(µ)0 ν1/3 (6.16)
and obtain
µ(ν, h) = 0.429(7)− 0.27(1)ν1/3, (6.17)
with χ2/d.o.f.= 1.1. This is in very good agreement with our balanced result. A plot
of the data and the fit is in Fig. 6.10.
6.5 Contact density
Another important quantity is the contact, which plays a prominent role in several
universal relations derived by Tan [59, 60, 61]. Through these relations the contact can
be interpreted in different ways, for instance as a measure of the local pair density [62], or
60
 0.08
 0.085
 0.09
 0.095
 0.1
 0.105
 0.11
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9
C/
ε F
2
ν1/3
C(ν = 0) = 0.1102(11)εF2
Figure 6.11: The contact density versus ν1/3 for the balanced data at the critical point. The
linear fit was performed for data at filling factors ν1/3 < 0.75 (filled circles).
Figure 6.12: Stereogram of the contact density C/ε2F versus filling factor and imbalance.
The surface corresponds to the quadratic fit (6.21). To view the image focus
your eyes on a point between the image and you (cross viewing).
61
a measure of the number of atoms with large momenta. The contact also features in the
expression for the total energy and in particular is responsible for making this expression
convergent, as the contact term cancels the divergent momentum distribution tail. The
adiabatic relation describes how the total energy changes with scattering length in terms
of the contact. In the presence of an external potential the contact appears in the virial
theorem and for a homogeneous system in the relation between pressure and energy.
We choose
C = m2g0Eint (6.18)
as the definition of the contact, where g0 denotes the physical coupling constant [62, 63].
The contact is related to the contact density C via
C =
∫
C(r)d3r, (6.19)
or for homogeneous systems simply C = CV . The dimensionless quantity C/ε2F can be
expressed through lattice quantities1 as
C/ε2F = (UEint)/(4L3ε2F ). (6.20)
Using only balanced data the best fit is C/ε2F = 0.1102(11)−0.033(2)ν1/3 with χ2/d.o.f. =
1.8. The data and the fit are shown in Fig. 6.11. In the presence of an imbalance we
expect the interaction energy and hence the contact density to decrease. A three di-
mensional fit of the data for C/ε2F to
C(ν, h) = C0 + C2h2 + (α(C)0 + α(C)2 h2)ν1/3 (6.21)
yields C0 = 0.1101(9), α(C)0 = −0.033(2), C2 = −0.15(16) and α(C)2 = −0.29(36), with
χ2/d.o.f.= 1.5. This is consistent with the balanced fit. The best fit value for the
parameter C2 is negative as required. Forcing α(C)2 = 0,
C(ν, h) = C0 + C2h2 + α(C)0 ν1/3, (6.22)
yields C0 = 0.1099(8), α(C)0 = −0.0322(15) and C2 = −0.01(2), with χ2/d.o.f.= 1.5. Fit-
ting to a constant function yields C(ν, h) = 0.1097(8)− 0.0320(14)ν1/3 with χ2/d.o.f.=
1.4. Figure 6.12 summarises our results. By a simple change of variables we can ob-
tain the contact in units of NkF from the contact density. Our most precise value for
the contact at zero imbalance (corresponding to C/ε2F = 0.1101(9) from fit (6.21)) is
C/NkF = 3.26(3).
1remember that we use m = 1/2 and hence in particular εF = k2F
62
 0
 0.05
 0.1
 0.15
 0  0.2  0.4  0.6  0.8  1
T c
/ε F
ν1/3
Figure 6.13: Linear continuum extrapolation of the balanced data for the critical tempera-
ture from this work [6] (red circles and red solid line) and the corresponding
extrapolation of the data from [26] (blue triangles and dashed blue line).
6.6 Comparison with the literature
Our final result for the critical temperature in physical units is Tc/εF = 0.171(5). This
value is significantly higher than the previous result from [26], where Tc/εF = 0.152(7).
In [64] the result of [26] was found to be in agreement with a continuous spacetime
DDMCmethod. As our method uses the same procedure as [26] to obtain the continuum
limit of the critical temperature we are able to directly compare our data at finite lattice
spacing, as shown in Fig. 6.13. From the plot we can see that the two data sets are
consistent with each other within statistical error. The difference in the continuum
result could be attributed to the different fitting range. In this work we only fit points
with ν1/3 < 0.75, while points at higher values of the filling factor were included in [26].
If the filling factor is too high second order lattice corrections become significant and
can distort the correct continuum result. In support of this, note that the three lowest
density points from [26] lie above the linear extrapolation.
The authors of [23] found an upper bound of Tc/εF / 0.15(1). They used an
auxiliary field Monte Carlo approach and extracted the critical temperature from the
finite-size scaling of the condensate fraction using the same procedure as in [26]. The
difference between our results might be attributed to the approximation made through
this fitting method, which was discussed in chapter 3. A less recent result by the
same group [24, 25] is Tc/εF = 0.23(2). Through extrapolating Monte Carlo results of
low-density neutron matter, the authors of [65] found a value of Tc/εF = 0.189(12) at
63
unitarity. Their value agrees with our result within errors. This work does not use the
finite-size scaling approach with systematic error, but determines Tc from the inflexion
point of the interpolated curve of the order parameter [66]. There are also results
obtained with the Restricted Path Integral Monte Carlo method [22], Tc/εF ≈ 0.245,
and an upper bound of Tc/εF < 0.14 obtained with a hybrid Monte Carlo method [21].
For completeness we also compare with several analytical predictions. A brief
overview over the methods was given in section 1.3. With the 1/N -expansion [15]
one obtains Tc/εF = 0.136 to leading order in 1/N . The beyond mean-field calculation
[20] yields Tc/εF = 0.225, however this result was matched with the less recent Monte
Carlo result from [24, 25], which has since then been updated to a lower value [23]. A
value of Tc/εF = 0.160 was obtained with the self-consistent approach [16]. Results
from the ²-expansion [17, 18] include Tc/εF = 0.249 using the expansion around d = 4
and several other values obtained from expansions around two and four dimensions with
different extrapolations. Together the results from these extrapolations range around
Tc/εF = 0.180(12), see [19] for a review. For comparison, the critical temperature in
the BEC limit is TBEC = 0.218εF .
The authors of [63] conjecture that the leading order change of the critical temper-
ature is linear in kF re, where re is the effective range of the potential, with a model
independent coefficient. Using the extrapolation of the balanced data our result for the
linear slope is ∆Tc/εF = −0.16(1), which is larger in magnitude than the value from
[26].
Our result for the energy per particle E/EFG = 0.440(15) shows excellent agreement
with the value E/EFG = 0.45(1) at Tc quoted in [23]. The value quoted in [26] is
E/NεF = 0.31(1), which roughly corresponds to E/EFG = 0.52(2). Our result for
the chemical potential µ/εF = 0.429(7) differs from µ/εF = 0.493(14) quoted in [26],
but is consistent with the value µ/εF = 0.43(1) quoted in [23]. Although the numerical
method used in [23] is different from our approach, our results for the energy per particle
and the chemical potential are in excellent agreement with each other. This supports
our claim that the discrepancy in the critical temperature between [23] and this work
is due to the systematic error in the finite-size scaling analysis of the order parameter
used in [23]. Analytical results for the chemical potential and the energy per particle
include µ/εF = 0.459 and E/NεF = 0.400 from [20], µ/εF = 0.394 and E/NεF = 0.304
from [16] and µ/εF = 0.18 and E/NεF = 0.212 from [18].
Since the chemical potential and the energy are expected to stay almost constant
at temperatures below Tc we also make a comparison to values from the literature
obtained at zero temperature. In the zero temperature limit the quantities µ/εF and
E/EFG are equal and Monte Carlo estimates range between approximately 0.40(1) and
0.44(1) [67, 68, 69, 70]. Our value for the chemical potential falls within this range, the
value for the total energy is slightly higher, which is consistent with the fact that the
energy must increase at finite temperature. These numerical estimates are consistent
64
 0.1
 0.12
 0.14
 0.16
 0.18
 0.2
 0.22
 0.24
 0.26
 0.28
T c
/ε F
Figure 6.14: Different predictions for the critical temperature at unitarity from the litera-
ture. Red circles denote numerical, blue squares analytical, and green triangles
experimental results. The solid line and the shaded area denote our most pre-
cise result and the error margin. From left to right the numerical results are:
the result from this work [6], the result from this work using only balanced
data [6], and the results from [26, 64, 23, 24, 65, 21, 22]. The analytical results
are from [20, 15, 16] and four results obtained with the ²-expansion [17, 18].
The experimental results are from [77, 78, 79]. For comparison the dashed line
denotes the critical temperature in the BEC limit.
with experiment [31, 71, 72].
Some predictions for the finite temperature contact density are also available, for
instance with a T -matrix approach [73] or by extracting the contact from the large
frequency tail of the shear viscosity [74]. The theoretical prediction [75] uses a nonper-
turbative virial expansion at high temperatures and compares different approximative
strong coupling theories at low temperatures. With the auxiliary field quantum Monte
Carlo approach the contact can be extracted from the high momentum tail of the den-
sity distribution [76]. Our result for the contact density at the critical point is in good
agreement with all these predictions. There are also results obtained at temperatures
much lower than Tc (see [62] and references therein).
Finally, we make a comparison with recent experimental studies of the homogeneous
unitary Fermi gas. A direct measurement of the critical temperature and the chemi-
cal potential of the uniform gas has been presented in [77]. Their experimental value
Tc/εF = 0.157(15) agrees well with our result. However, the value of the chemical po-
tential at the critical point µ/εF = 0.49(2) differs from our value. Another experimental
determination of the critical temperature and thermodynamic functions, including the
65
energy and the chemical potential, is described in [78]. Their values Tc/εF = 0.17(1)
and µ/εF = 0.43(1) at Tc show excellent agreement with our results. Their result for
the energy per particle E/NεF = 0.34(2) at Tc is higher than our value. In another
experimental work [80] an estimate for the critical temperature at zero imbalance is
extrapolated from data at higher values of imbalance. In the most recent experimental
study of the superfluid transition of a homogeneous unitary Fermi gas [79] the critical
temperature is determined from the sudden rise of the specific heat at the critical point,
yielding Tc/εF = 0.167(13), which agrees very well with our result. Their value for the
chemical potential at Tc is µ/εF = 0.42(1), which is also consistent with our result. An
experimental value for the contact density of the homogeneous unitary Fermi gas at
zero temperature, C/ε2F = 0.1184(64), is presented in [81]. In Fig. 6.14 we compare the
different values of the critical temperature found in the literature with our results.
66
Chapter 7
Observables beyond the critical point
So far we have been limited to observables calculated at the critical point. Now we want
to study how the chemical potential, the energy per particle and the contact density
change with temperature. Although the lattice temperature T = 1/β is a tunable
parameter, fixing the physical temperature for different filling factors is a non-trivial
task. The critical point is special in this respect, as it has a unique physical property
that allows us to distinguish it from other points in the phase diagram. By looking at the
behaviour of the order parameter we can identify the critical lattice temperature for any
given value of the lattice chemical potential µ, and hence create a linear extrapolation,
for instance like the one presented in Fig. 6.2. Hence it is convenient to use the lattice
critical temperature Tc(µ) as a benchmark to fix the temperature scale.
When working on the lattice, in addition to the physical length scale determined
by the interparticle distance we have another (artificial) length scale given by the lat-
tice spacing b. The lattice spacing is responsible for the non-zero slope of the linear
extrapolations of physical observables and is eventually removed through taking the
continuum limit. Remember from chapter 2 that we set the lattice spacing via ν = nb3.
For instance we could have used an experimental value for the density n of an atomic
gas at the critical point to compute b at any given filling factor in physical units. In
practice, due to the universality of the unitary Fermi gas, a matching with experimental
values is not necessary.
Now let us vary the temperature. The simplest approach is to hold the lattice
chemical potential fixed and to vary T = 1/β. We now set the lattice spacing such that
it is independent of T ,
b(µ, T ) = b(µ, Tc) =
(ν(µ, Tc)
n(µ, Tc)
)1/3
. (7.1)
This can also be understood as a temperature-independent renormalisation condition. If
we fix the lattice temperature ratio r = T (µ)/Tc(µ) for each value of the lattice chemical
potential, we will move along a line of constant temperature T = rTc analogous to the
linear extrapolations at Tc discussed in the previous chapter.
For coarse lattices this scheme will break down as the presence of the lattice spacing
will change the relation between ν and T . We need to ensure that for the temperature
67
 0
 2
 4
 6
 8
 10
 0  1  2  3  4  5  6  7  8  9  10
ν(µ
,
T)/
ν(µ
,
T c
)
T/Tc
µ=0.4
µ=0.5
µ=0.7
Figure 7.1: The normalised filling factor ν(µ, T )/ν(µ, Tc) versus the temperature ratio T/Tc
for different values of the lattice chemical potential µ.
range considered in this analysis the lattice artefacts are negligible. Note that the filling
factor normalised by its value at the critical point must be independent of the lattice
chemical potential. Figure 7.1 shows the normalised filling factor versus T/Tc for several
values of µ. It is clearly visible that we can access temperature ratios T/Tc ≤ 4 provided
that the chemical potential is small enough. For instance the chemical potential value
µ = 0.7 in lattice units is too large to access the temperature ratio T/Tc = 4, but small
enough for T/Tc ≤ 3. The deviation becomes more evident for very large temperature
ratios. In general, high temperatures are hard to access as they require us to go to very
low chemical potentials and consequently very large lattice sizes.
In the following sections we will study the temperature dependence of the chemical
potential, the energy per particle and the contact density of the balanced unitary Fermi
gas. Data was taken below the critical temperature at T/Tc = 0.7 and for six different
temperature ratios above the critical temperature, up to T/Tc = 4. This analysis is
still work in progress and the results presented in the following are preliminary. In
particular, simulations at the lowest filling factor that was used for the analysis at
Tc are still ongoing. At high temperature and low filling finite-size effects are more
pronounced and force us to use large lattice sizes which require long simulation times.
68
 0.1
 0.15
 0.2
 0.25
 0.3
 0.35
 0.4
 0.45
 0.5
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
µ/
ε F
ν1/3
T/Tc = 4.0T/Tc = 3.5T/Tc = 3.0T/Tc = 2.5T/Tc = 2.0T/Tc = 1.5T/Tc = 1.0T/Tc = 0.7
Figure 7.2: The chemical potential versus filling factor for different temperatures. The con-
tinuum limit corresponds to ν → 0. The solid line is the continuum extrapolation
at Tc.
T/Tc linear fit for µ/εF χ2/d.o.f.
4.0 0.20(2)− 0.08(3)ν1/3 0.76
3.5 0.25(2)− 0.14(4)ν1/3 0.16
3.0 0.25(1)− 0.10(2)ν1/3 1.37
2.5 0.29(1)− 0.15(2)ν1/3 3.44
2.0 0.35(1)− 0.20(2)ν1/3 1.67
1.5 0.41(2)− 0.27(3)ν1/3 1.12
0.7 0.38(1)− 0.20(2)ν1/3 0.80
Table 7.1: Linear fits of the numerical data for the chemical potential µ/εF at different
temperatures.
7.1 The equation of state for the density
Figure 7.2 displays the numerical data for the chemical potential at different tempera-
tures. For lucidity the linear fits at T 6= Tc are not shown, but a list of all fit results
is given in Table 7.1. It is clearly visible from the plot that the chemical potential
decreases with increasing temperature. Figure 7.3 shows the continuum limit of the
chemical potential as a function of the temperature.
For comparison with the literature we will now rewrite the equation of state in terms
of the density as a function of βµ. At unitarity the dimensionless density nλ3B, where
λB =
√
2pi/mT is the thermal de Broglie wavelength, must be a universal function of
69
 0.1
 0.2
 0.3
 0.4
 0.5
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
µ/
ε F
T/εF
Figure 7.3: The chemical potential in the continuum limit versus the temperature. The red
symbol indicates the value at the critical point.
 1
 2
 3
 4
-2 -1  0  1  2  3  4
f n(
βµ
)/f n
(0)
(βµ
)
βµ
Figure 7.4: The density fn(βµ)/f (0)n (βµ) normalised in terms of the density of the ideal
gas. We compare our results (red circles) with preliminary results obtained
with bold diagrammatic Monte Carlo [82] (blue triangles) and the third order
virial expansion [83] (solid line).
70
βµ,
nλ3B = fn(βµ). (7.2)
Using the definition of the Fermi energy εF = (3pi2n)2/3/(2m) we can rewrite
fn(βµ) = nλ3B =
8
3√pi (T/εF )
−3/2. (7.3)
For an ideal two-component Fermi gas this universal function can be evaluated explicitly
from the Fermi-Dirac distribution,
f (0)n (βµ) = n(0)λ3B = 2
∫ d3p
(2pi)3
1
1 + eβ( p22m−µ)
( 2pi
mT
)3/2
= −2Li3/2(−eβµ), (7.4)
where Lin(z) =
∑∞
k=1(zk/kn) is the polylogarithm. Figure 7.4 shows our results for the
density normalised in terms of the density of the ideal gas. We also compare with the
unpublished preliminary results obtained using bold diagrammatic Monte Carlo1 [82]
and with a theoretical prediction valid at high temperatures obtained from the virial
expansion (the power series expansion of thermodynamic quantities in terms of the
fugacity eβµ). For theoretical background on the virial expansion of the unitary Fermi
gas and the second virial coefficient see for instance [85] and references therein. The
third virial coefficient for the unitary Fermi gas was calculated in [83]. We can see good
agreement between our results and the results from [82]. At very high temperatures
lattice artefacts become more pronounced and the results start to deviate. The results
obtained with the auxiliary field quantum Monte Carlo method [23] (not shown in the
plot) also agree with our predictions.
7.2 The equation of state for the pressure
The numerical results for the energy per particle for different temperatures are shown
in Fig. 7.5. Again only the linear extrapolation at Tc is shown and the fit results are
summarised in Table 7.2. As expected the energy per particle increases with increasing
temperature. A plot of E/EFG in the continuum limit versus the temperature is shown
in Fig. 7.6.
Using universal thermodynamic relations we can express the pressure P of the gas
in terms of the temperature, the chemical potential and the energy density. On dimen-
sional grounds we obtain
Pλ3B = TfP (βµ), (7.5)
where fP (βµ) is a universal function. The pressure can be related to the energy density
1A more recent version of these results and a description of the method can be found in [84].
71
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
E/
E F
G
ν1/3
T/Tc = 4.0T/Tc = 3.5T/Tc = 3.0T/Tc = 2.5T/Tc = 2.0T/Tc = 1.5T/Tc = 1.0T/Tc = 0.7
Figure 7.5: The energy per particle versus filling factor for different temperatures. The
continuum limit corresponds to ν → 0. The solid line is the continuum extrap-
olation at Tc.
T/Tc linear fit for E/EFG χ2/d.o.f.
4.0 1.08(5)− 0.89(7)ν1/3 4.45
3.5 1.03(5)− 0.82(8)ν1/3 0.48
3.0 1.06(4)− 0.92(6)ν1/3 1.21
2.5 0.96(3)− 0.79(5)ν1/3 1.76
2.0 0.86(4)− 0.68(6)ν1/3 2.39
1.5 0.79(4)− 0.63(8)ν1/3 1.53
0.7 0.39(4)− 0.18(7)ν1/3 4.58
Table 7.2: Linear fits of the numerical data for the energy per particle E/EFG at different
temperatures.
of the gas [47] via P = (2/3)E/V . From this we get
fP (βµ) = Pλ
3
B
T =
2Eλ3B
3V T =
16
15√pi (E/EFG)(T/εF )
−5/2. (7.6)
Again we would like to normalise the pressure in terms of the pressure of the ideal
two-component Fermi gas,
f (0)P (βµ) =
2E(0)λ3B
3V T =
4
3T
∫ d3p
(2pi)3
p2/2m
1 + eβ( p22m−µ)
( 2pi
mT
)3/2
= −2Li5/2(−eβµ). (7.7)
Figure 7.7 shows our results for the equation of state for the pressure. For comparison
72
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
E/
E F
G
T/εF
Figure 7.6: The energy per particle in the continuum limit versus the temperature. The red
symbol indicates the value at the critical point.
 0
 1
 2
 3
-2 -1  0  1  2  3  4
f P
(βµ
)/f P
(0)
(βµ
)
βµ
Figure 7.7: The pressure fP (βµ)/f (0)P (βµ) normalised in terms of the pressure of the ideal
gas. We compare our results (red circles) with the third order virial expansion
[83] (solid line).
73
 0.09
 0.1
 0.11
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
C/
ε F
2
ν1/3
T/Tc = 4.0T/Tc = 3.5T/Tc = 3.0T/Tc = 2.5T/Tc = 2.0T/Tc = 1.5T/Tc = 1.0T/Tc = 0.7
Figure 7.8: The contact density versus filling factor for different temperatures. The contin-
uum limit corresponds to ν → 0. The solid line is the continuum extrapolation
at Tc.
we also show the high temperature prediction from the virial expansion. The agreement
here is less good than in the case of the equation of state for the density. We expect an
improvement of the data when results at the lowest filling factor become available.
7.3 The temperature dependence of the contact
Figure 7.8 shows our data for the contact density at different temperatures and the
linear fit at Tc. Like the energy per particle the contact apprears to increase with
increasing temperature. The fit results are listed in Table 7.3. Figures 7.9 and 7.10
show the contact density versus the temperature and versus βµ in the continuum limit.
We also compare our results to the literature. There is a discrepancy to the preliminary
results from [82] where the contact was found to be increasing with increasing βµ. On
the other hand we see excellent agreement with the auxiliary field quantum Monte Carlo
results [76] that show the same trend of the contact density with respect to βµ as our
results. We also agree with the T -matrix calculation in [73]. Our result at the lowest
temperature value T/Tc = 0.7 shows excellent agreement with the zero-temperature
experimental result [81].
74
 0.1
 0.11
 0.12
 0.13
 0.14
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8
C/
ε F
2
T/εF
Figure 7.9: The contact density in the continuum limit versus the temperature. The red
symbol indicates the value at the critical point.
 0.05
 0.06
 0.07
 0.08
 0.09
 0.1
 0.11
 0.12
 0.13
 0.14
-4 -3 -2 -1  0  1  2  3  4
C/
ε F
2
βµ
Figure 7.10: The contact density in the continuum limit versus βµ. Red circles are the
results from this work, blue triangles the preliminary results from [82] and
the lines indicate the experimental zero-temperature result [81] with the error
margin.
75
T/Tc linear fit for C/ε2F χ2/d.o.f.
4.0 0.134(4)− 0.067(7)ν1/3 0.01
3.5 0.131(4)− 0.062(7)ν1/3 0.72
3.0 0.127(3)− 0.056(5)ν1/3 0.66
2.5 0.124(2)− 0.051(3)ν1/3 1.43
2.0 0.115(2)− 0.038(3)ν1/3 1.79
1.5 0.107(2)− 0.025(3)ν1/3 1.16
0.7 0.115(2)− 0.040(4)ν1/3 0.43
Table 7.3: Linear fits of the numerical data for the contact density C/ε2F at different tem-
peratures.
76
Chapter 8
Simulation of the Boltzmann equation
The previous chapters were devoted to the homogeneous unitary Fermi gas at or near
the critical point. Now we will move to higher temperatures, when a semiclassical de-
scription is applicable. We shall also include an external potential and consider a range
of large but finite values for the scattering length a. Our goal is to access dynamical
properties, in particular related to spin transport, using a numerical simulation of the
Boltzmann equation.
8.1 Distributions in a trap
Again we consider a system of two-component fermions with equal mass m, labelled by
the spin index s = {↑, ↓}. As we have seen in section 1.2, fermions of opposite spin can
interact via s-wave scattering and the cross section is given by (1.13),
σ = 4pia
2
1 + k2a2 , (8.1)
where 2k = |p↑ − p↓| = |prel| is the relative momentum of the two atoms.
We assume that the system is in the normal phase and that the temperature is
sufficiently high, so that the two spin distributions can be described semiclassically
in terms of functions fs(r,p, t) [3, 13]. The distribution function fs(r,p, t) gives the
probability density to find a particle with spin s at the phase space point (r,p) at
the time t. Within the local density approximation the equilibrium distribution for
fermions is given by the Fermi-Dirac distribution
f (FD)s (r,p) =
1
e(p2/2m+V (r)−µs)/T + 1 , (8.2)
where T is the temperature and µs the chemical potential defined by the normalisation
condition
Ns = N/2 =
∫
d3rns(r, t) =
∫
d3r
∫ d3p
(2pi)3fs(r,p, t). (8.3)
The atom number for each species is assumed to be equal, N↑ = N↓ = N/2, and will be
77
held fixed during the simulation. The trapping potential is assumed to be harmonic,
V (r) = 12m(ω
2
xx2 + ω2yy2 + ω2zz2), (8.4)
with the three trapping frequencies ωx, ωy and ωz. In the following we will only consider
the isotropic case ωx = ωy = ωz or a cigar shaped trap with ωz < ωx = ωy.
It will prove useful to derive the corresponding energy distribution by performing a
change of variables,
f (FD)s (E) =
∫
d3r d
3p
(2pi)3f
(FD)
s (r,p)δ(E − E(r,p)) (8.5)
=
∫
E(r,p)=E
d3r d
3p
(2pi)3
f (FD)s (r,p)√∑
i (∂E/∂ri)2 +
∑
i (∂E/∂pi)2
, (8.6)
where E(r,p) = p2/2m+V (r) is the total energy. Performing the integral in momentum
space we then obtain
f (FD)s (E) =
1
e(E−µs)/T + 1
m
pi2
∫
E−V (r)≥0
d3r E − V (r)√2(E − V (r))/m+m2∑i ω4i r2i
. (8.7)
On dimensional grounds and given the integration condition E − V (r) ≥ 0 it is now
easy to see that the integral must be proportional to E2. From normalisation we obtain
f (FD)s (E) =
E2/2ω30
e(E−µs)/T + 1 , (8.8)
where ω0 = (ωxωyωz)1/3 is the geometric average over the three trapping frequencies.
The chemical potential at zero temperature defines the Fermi energy,
ε˜F = µs(T = 0), (8.9)
which in the presence of interactions differs from the Fermi energy εF of the free non-
interacting gas [3]. An explicit expression for ε˜F can be obtained from the normalisation
of the energy distribution at zero temperature,
Ns = limT→0
1
2ω30
∫ ∞
0
E2dE
e(E−ε˜F )/T + 1 =
1
2ω30
∫ ε˜F
0
E2dE = ε˜
3
F
6ω30
, (8.10)
which yields
T˜F = ε˜F = (6Ns)1/3ω0 = (3N)1/3ω0. (8.11)
For better comparison with the experiment we will sometimes also use εF as a unit of
energy. Since in the presence of an external potential the gas is no longer homogeneous,
78
we need to redefine εF in terms of a local density. We take
TF = εF = 12m(3pi
2n(0))2/3, (8.12)
where n(0) is the total atomic density n(0) = n↑(0) + n↓(0) in the trap centre. The
local density of the trapped gas is temperature-dependent and consequently εF and
kF =
√2mεF defined in terms of n(0) will change if the temperature of the system
changes. The zero-temperature Fermi energy ε˜F is defined in terms of only the total
particle number and the trap geometry and stays invariant if we hold these quantities
fixed. Note that in contrast to the dilute homogeneous gas which can be described
entirely in terms of the density, the temperature and the scattering length, the trapping
frequencies introduce three additional scales into the problem.
The zero-temperature Fermi energy ε˜F can be used to determine the typical scales
of the cloud,
ε˜F = k˜
2
F
2m =
1
2mω
2
xR2x =
1
2mω
2
yR2y =
1
2mω
2
zR2z, (8.13)
where k˜F is the zero-temperature Fermi momentum and Rx, Ry and Rz are the Thomas-
Fermi radii in the three spatial directions. The Fermi momentum and the Thomas-
Fermi radii give the width of the zero-temperature momentum and density distributions
respectively. They are therefore useful quantities to describe the extent of the fermionic
cloud in momentum and coordinate space. To see the correspondence between these
quantities and the widths of the distributions let us look, for instance, at the Fermi-
Dirac density distribution in the zero-temperature limit,
lim
T→0
n(FD)s (r) = limT→0
∫ d3p
(2pi)3
1
e(p2/2m+V (r)−µs)/T + 1 =
∫
p2/2m+V (r)≤ε˜F
d3p
(2pi)3 . (8.14)
Using rotational symmetry we obtain
lim
T→0
n(FD)s (r) =
∫ √2m(ε˜F−V (r))
0
4pip2dp
(2pi)3 =
1
6pi2 (2m(ε˜F − V (r)))
3/2 (8.15)
and inserting the definition of the Thomas-Fermi radii from (8.13) we arrive at
lim
T→0
n(FD)s (r) =
8Ns
pi2RxRyRz
(
1− x
2
R2x
− y
2
R2y
− z
2
R2z
)3/2
. (8.16)
Analogously
lim
T→0
n(FD)s (p) = limT→0
∫
d3rf (FD)s (r,p) =
8Ns
k˜3F
(
1− p
2
k˜2F
)3/2
. (8.17)
79
Note that unlike the density distribution the momentum distribution is isotropic in the
three directions, even if the trap frequencies are not equal.
For sufficiently high temperatures the exponential e−(p2/2m+V (r)−µs)/T becomes small
such that we can expand the Fermi-Dirac distribution in this quantity and obtain
f (FD)s (r,p) =
e−(p2/2m+V (r)−µs)/T
1 + e−(p2/2m+V (r)−µs)/T → e
−(p2/2m+V (r)−µs)/T , (8.18)
which is the classical Maxwell-Boltzmann distribution. Using the normalisation condi-
tion (8.3) we can rewrite it explicitly in terms of the particle number rather than the
chemical potential,
f (MB)s (r,p) = Ns
ω30
T 3 e
−(p2/2m+V (r))/T . (8.19)
The corresponding energy distribution is
f (MB)s (E) = Ns
E2
2T 3 e
−E/T . (8.20)
The classical Maxwell-Boltzmann distribution is applicable if the de Broglie wavelength
is much smaller than the average interparticle distance, which corresponds to T/T˜F & 1.
8.2 The Boltzmann equation
Within the scope of kinetic gas theory a many-body system is treated as a collection
of discrete particles moving randomly [49, 86]. If the system is sufficiently dilute such
that the size of the particles is negligible compared to the interparticle separation,
the particles can be viewed as structureless point-like objects. They will propagate
freely most of the time, interparticle interactions being limited to rare elastic two-body
collisions. If the de Broglie wavelength is sufficiently small the free propagation follows
the classical equations of motion.
We are interested in the non-equilibrium properties of the system and hence want
to look at the time-evolution of the distribution function fs at the phase space point
(r,p). For a closed system the distribution function can only change due to collisional
processes. Hence its total time derivative can be written as
d
dtfs =
∂
∂tfs + p˙ · ∇pfs + r˙ · ∇rfs = −I, (8.21)
where the quantity I represents the effect of two-body collisions and in general de-
pends on the two-body distribution function f (2)ss (rs,ps, rs,ps, t). The analogous time-
evolution equation for f (2) will in turn involve the three-body distribution f (3) and so
on [86]. In order to obtain a closed form for the time-evolution equation we need to
make an additional assumption. In a dilute gas of point-like particles with short-range
80
interactions collisions are rare and hence it is justified to assume that before a collision
the two particles are uncorrelated,
f (2)ss (rs,ps, rs,ps, t) = fs(rs,ps, t)fs(rs,ps, t). (8.22)
This is called the molecular chaos assumption or Stoßzahlansatz. After the collision
of course, the two particles become strongly correlated due to momentum and energy
conservation. But since they are unlikely to encounter each other again before scattering
on other particles has erased the memory of their prior collision the molecular chaos
assumption remains valid. More precisely, it can be shown that the probability for such
a repeated scattering event is exponentially small with time [86].
Now let us derive an explicit form for the collisional contribution I taking into
account all assumptions discussed above [49]. The change of the distribution function
fs at the phase space point (r,p) can be decomposed into a gain and a loss term, the
gain term representing atoms of spin s with final state (r,p) after the collision and the
loss term representing the contribution of collisions from the initial spin s state (r,p)
into other final states (r,p′). Now consider a specific collision between a spin s atom in
the initial state (r,ps) and an atom of opposite spin s in the initial state (r,ps), which
brings the two atoms into the final states (r,p′s) and (r,p′s) respectively. Collisions are
assumed to be local, meaning that the two atoms are at the same position r during
the collision. The total number of such collisions with fixed initial and final momenta
divided by the phase space volume is given by the initial particle density fs(r,ps) times
the probability for each spin s atom to exhibit such a collision. The latter in turn
equals the product of the number of possible scattering partners fs(r,ps)d3ps with the
probability to end up in the required final state, |vs − vs|dσ = |vs − vs| dσdΩdΩ, where
Ω is the solid angle between the incoming and outgoing relative velocities as discussed
in section 1.2. The change of fs due to all such collisions becomes an integral over all
possible initial momenta of the collision partners ps and over all possible angles Ω. For
“classical fermions” this collision integral then takes the form
Iclass[fs, fs] =
∫ d3ps
(2pi)3
∫
dΩ dσdΩ
|ps − ps|
m [fsfs − f
′
sf ′s], (8.23)
where the primed variables refer to quantities after the collision. Since fermions obey
the Fermi-Dirac statistics a particle can only scatter into a previously unoccupied quan-
tum state. This reduces the scattering probability by a so-called Pauli blocking term
proportional to (1− f ′s)(1− f ′s). Taking this into account the collision integral reads
I[fs, fs] =
∫ d3ps
(2pi)3
∫
dΩ dσdΩ
|ps − ps|
m [fsfs(1−f
′
s)(1−f ′s)−f ′sf ′s(1−fs)(1−fs)]. (8.24)
If we now insert the explicit form for the collision integral and the classical equations
81
of motion into (8.21) we arrive at the Boltzmann equation which describes the time-
evolution of the distribution function fs(r,p, t) for each spin species,
∂tfs + (p/m) · ∇rfs −∇rV · ∇pfs = −I[fs, fs]. (8.25)
The left-hand side represents the propagation of the atoms in the potential while the
right-hand side stands for the collision integral given by Eq. (8.24) for fermions or by
Eq. (8.23) for classical particles.
8.3 Numerical setup
The results presented in the following were published in [8].
Our numerical setup is based to the one described in [87]. We do not attempt to
calculate the collision integral (8.24) directly, as this is numerically intractable. Instead
we start from the semiclassical picture of a dilute gas of point-like particles propagating
in space and scattering on each other via rare elastic two-body collisions. We introduce
a discrete time step ∆t, such that during each time step the atoms propagate collision-
free following their classical trajectories. At the end of each time step collisions between
the atoms are evaluated. The point-like particle picture is a discrete approximation of
the continuous distribution function fs(r,p, t) through δ-functions, which allows us to
handle the propagation step efficiently. In order for this approximation to be accurate
we will represent each fermion by several test particles [87, 88, 89]. The higher the ratio
N˜/N of test particles to atoms, the more precisely the continuous distribution function
will be approximated,
fs(r,p, t) → NsN˜s
N˜s∑
i=1
(2pi)3δ(r− ri(t))δ(p− pi(t)). (8.26)
Physical observables need to be rescaled, for instance the test particle cross section
becomes σ˜ = σ(N/N˜). A generic thermal expectation value of a single-particle observ-
able X(r,p) can be easily calculated within the test particle picture, as the integration
reduces to a sum over all test particles,
〈X〉 =
∑
s
1
Ns
∫
d3r d
3p
(2pi)3fs(r,p, t)X(r,p) =
1
N˜
N˜∑
i=1
X(ri,pi). (8.27)
The trajectories of the test particles are the same as the trajectories of the atoms and are
82
given by the stationary solution of the Boltzmann equation with a harmonic potential,
ri(tn+1) = ri(tn) cos(ωi∆t) + (pi(tn)/mωi) sin(ωi∆t), (8.28)
pi(tn+1) = pi(tn) cos(ωi∆t)− ri(tn)mωi sin(ωi∆t). (8.29)
Note that since the time step is fixed the trigonometric functions only need to be
evaluated once during the entire simulation, so that using the exact solution is more
efficient than using the Verlet algorithm, as for instance in Refs. [87, 88], in which case
the accelerations need to be recalculated for every time step. The Verlet algorithm
is more general as it is applicable for any potential. Since in this work we will only
consider harmonic potentials, we will use the exact solution which is both more precise
and more efficient to calculate.
We evaluate collisions in the same way as described in [87]. First we test whether a
pair of test particles reaches the point of closest approach during the present time step.
This condition is important to prevent particles from attempting to collide with each
other repeatedly for several consecutive time steps. If the closest approach condition
is true we check if the minimal distance dmin of the test particles fulfils the classical
condition for scattering: pid2min < σ˜. If this condition is also satisfied we propose a
collision at the time of closest approach. However due to Pauli statistics, even if the
classical conditions for scattering are fulfilled, a collision can only take place if the new
state of the particles was previously unoccupied. To take this into account we calculate
the quantum mechanical scattering probability given by the Pauli term (1 − f ′s)(1 −
f ′s) in the collision integral (8.24) and accept or reject the collision according to this
probability. Clearly the point-like particle picture is unsuitable for the calculation of
this probability. To return to a continuous distribution we therefore have to smear
out the δ-functions representing the test particles, e.g. by Gaussians in position and
momentum space:
δ(p− pi)δ(r− ri) −→ e
−(p−pi)2/w2p
(√piwp)3
e−(x−xi)2/w2x√piwx
e−(y−yi)2/w2y√piwy
e−(z−zi)2/w2z√piwz . (8.30)
The widths of these Gaussians, wx, wy, wz and wp, need to be tuned so that on the one
hand fluctuations due to the discrete nature of the test particle picture are smoothed
out, but on the other hand the physical structure of the distribution function fs remains
preserved [87]. The statistical fluctuations are of the order of (N˜w3pwxwywz/N)−1/2, so
that the first condition is equivalent to
wpwr À (N/N˜)1/3, (8.31)
where we introduced wr = (wxwywz)1/3, the geometric average of the spatial widths.
The second condition implies wi ¿ Ri and wp ¿ k˜F , where the Fermi momentum
83
and the Thomas-Fermi radii were defined in (8.13). Additionally at low temperature
it is crucial to resolve the rapid change of the distribution function around the Fermi
surface, such that
wp ¿ k˜F (T/T˜F ) and wi ¿ Ri(T/T˜F ). (8.32)
Note that the smearing width in momentum space is isotropic, while in position space
the smearing width can be different depending on the spatial direction, if the corre-
sponding trap frequencies are unequal. This is a direct consequence of the isotropy of
the momentum space distribution (8.17) that was discussed in section 8.1. Since the
Thomas-Fermi radii are inversely proportional to the corresponding trap frequencies
it is sensible to choose wi = wrω0/ωi for the spatial widths. Furthermore Eq. (8.32)
together with the definition (8.13) imply wp = mwr. Hence all four smearing widths
can be reduced to only one free parameter. At very low temperatures the margin given
by the conditions (8.31) and (8.32) becomes so narrow that it is impossible to find
smearing widths satisfying both conditions, without having to significantly increase the
number of test particles. This limits the applicability of this setup to temperatures
above approximately 0.2T˜F for N˜/N = 10. Moreover, at very low temperatures the
system undergoes a phase transition into a superfluid state. This algorithm does not
include the relevant degrees of freedom of the superfluid Fermi gas and is only applicable
in the normal phase.
The main numerical challenge is to efficiently evaluate collisions between the test
particles. The total number of pairs of opposite spin is N˜2/4, which is an unfavourable
scaling given that we want to use large particle numbers and a high test particle to
particle ratio. In this work we develop a more efficient method than to check all possible
test particle pairs. The key observation is that since the cross section is decreasing with
increasing relative momentum it can never exceed σ˜max = 4pia2N/N˜ and consequently
the maximal distance two colliding test particles can have is dmax = 2a
√
N/N˜ . Having
this in mind we superpose a cubic grid with cell size dmax on the continuous space. The
grid has finite extent which can be set by demanding that a certain proportion, for
instance at least 95%, of the test particles are within the grid. For a cigar shaped trap
the grid must have larger extent in the axial direction. At the end of each time step
we move systematically through all grid cells starting in one of the corners and note all
possible collision partners that fulfill the classical scattering conditions. To make sure
that we do not miss collisions due to boundary effects, for each particle we check not
only all opposite spin particles in the same grid cell, but also in all neighbouring cells
(the ones sharing a face, an edge or a vertex with the given cell). This ensures that all
particles in a sphere of radius dmax around a given particle are definitely accounted for.
This makes a total of 33 = 27 cells for each particle, however to avoid double counting
we only need to evaluate cells in the positive direction.
A small systematic error source remains with this setup. If the relative velocity of
84
two particles is large they can be in non-neighbouring cells at the beginning and at
the end of a time step, although in the course of the time step they come within each
other’s allowed collision range. Such a possible collision will then not be accounted
for. However this systematic error can be minimised by choosing the time step to be
sufficiently small and also by choosing the cell size to be larger than dmax. Also note
that for large relative velocities the cross section is small and collisions between very
fast particles are rare events.
After having searched the entire grid for classically allowed collision pairs we proceed
to choose which collisions will indeed take place. To do so we consecutively select
random pairs from the list. We then propagate both particles to the point of their
closest approach, let them scatter (the exact setup for determining the new positions
and momenta after scattering is described below) and then propagate them back to the
original time. To account for quantum statistics we then calculate the Pauli blocking
factors using the new distributions f ′s = fs(r′,p′). The continuous distributions are
obtained by replacing the δ-functions in (8.26) according to (8.30),
fs(r′,p′) = NsN˜s
N˜s∑
i=1
(2pi)3 e
−(p′−pi)2/w2p
(√piwp)3
e−(x′−xi)2/w2x√piwx
e−(y′−yi)2/w2y√piwy
e−(z′−zi)2/w2z√piwz . (8.33)
If we obtain a value f ′s > 1 from the summation we set f ′s = 1. The probability that
the collision is accepted is then given by (1 − f ′s)(1 − f ′s). Regardless if a collision is
accepted or rejected neither of the particles concerned is allowed to collide again with
another particle during the present time step. If the collision is accepted we keep the
new positions and momenta. If the collision is rejected we return to the values before
the collision. This procedure is repeated until all possible pairs have been evaluated.
To calculate the new coordinates and momenta of the two particles after a collision
we use the following procedure. As collisions are local, disregard for the moment the ex-
ternal potential and go to the centre of mass frame of the two particles. Motivated by the
analogue of classical scattering we wish to conserve the total momentum, the total en-
ergy and the total angular momentum of the system during the collision. Conservation
of the total momentum p = ps + ps and the total energy E = Ekin = (p2s + p2s)/2m to-
gether imply the conservation of the modulus of the relative momentum prel = |ps−ps|,
since p2 + p2rel = 2(p2s + p2s) = 4mE. The direction of the relative momentum vector
can change during the collision. Conservation of the angular momentum L = rrel×prel
implies that the relative momentum vector can only be rotated in the plane spanned
by rrel and prel, or in other words prel must be rotated around the axis defined by the
vector L. The angle of this rotation is the only degree of freedom of the collision and is
determined at random. To conserve the modulus of the angular momentum, the relative
position vector rrel must be rotated by the same angle as prel. The relative distance of
the particles |rrel| still remains unchanged. It is easy to see that the angular momen-
85
tum and the relative position vector cannot be cannot be simultaneously conserved. We
chose a setup in which momentum, energy and angular momentum are conserved at the
cost of a rotation of the relative position vector of the particles. From the new values
p′rel and r′rel the new values of the momenta and positions of the individual particles can
be recovered using total momentum conservation and the centre of mass coordinates
respectively.
So far we have ignored the external potential, which is justified since the collisions
are local and hence the presence of an external potential should not play a role for
the calculation of the new coordinates and momenta. In practice however, in our
setup the two colliding particles are not exactly at the same position and their relative
position changes after a successful collision. Thus potential energy is not conserved and
energy conservation is not exact during a collision if the potential is non-isotropic. On
average of course, energy conservation still holds since in a many-particle system these
small effects will cancel each other out. It is possible to preserve energy conservation
exactly by employing a different setup, for instance as in [87]. In this reference, the
relative position stays fixed during a collision and the relative momentum vector is
rotated by a random angle in space (such a rotation has two degrees of freedom). As a
direct consequence of the unrestricted rotation this setup violates angular momentum
conservation. On the other hand it is arguable that angular momentum is not conserved
during the propagation step either if the external potential is non-isotropic. As in this
work we will always consider either isotropic or axially symmetric potentials, either the
total angular momentum or its axial component Lz are conserved.
To study the impact of the collisional setup on the outcome of the simulation we
have implemented the setup from [87] in addition to our setup described above, and
confirmed that they generate the same results for both equilibrium and non-equilibrium
systems. We compared the equilibrium collision rates as well as the time evolutions of
the centre of mass coordinates in response to a perturbation. The differences between
the two setups were found to be smaller than the statistical fluctuations. It is therefore
well-grounded to say that the two setups are equivalent on a macroscopical scale and we
have the freedom to choose whichever option appears more convenient. We adopt the
angular momentum conserving setup, not only because of this property, but also since
the setup from [87] has a small technical disadvantage. Since the new direction of the
relative momentum vector is chosen uniformly on a sphere, in many cases the particles
are found to approach each other again after a successful collision. This implies that
in the following time step they are likely to undergo another collision. To avoid this
overcounting of collisions one needs to implement an additional routine that prevents
particles from colliding with each other repeatedly within short time intervals. In other
words the molecular chaos assumption does not hold with the setup from [87] and needs
to be enforced artificially. In our setup repeated scattering is so rare that this effect
can be neglected.
86
Since the picture of colliding point-like particles with well-defined positions and
momenta is a classical interpretation of a quantum mechanical scattering process, it
is unsurprising that there is some ambiguity in the implementation of the collisional
setup. In the semiclassical particle picture each collision has 12 degrees of freedom:
three position and three momentum components for each of the two particles, or equiv-
alently three components of the total and relative positions and momenta. In quantum
mechanics we consider wave packets rather than particles and concepts like particle po-
sition or momentum are not well-defined. Instead the system is described by operators
which, depending on the symmetries of the system, commute with the Hamiltonian
and define quantum numbers corresponding to conserved quantities. The concept of a
trajectory does not exist, as a particle is not localised in phase space but rather smeared
out over a certain phase space volume in accordance with Heisenberg’s uncertainty prin-
ciple. It is therefore not necessary to preserve for instance the exact positions of the
two atoms during a collision. In the presence of an axially symmetric external potential
for instance, the only conserved quantities are the total energy and Lz.
From this quantum mechanical picture we can develop the following general colli-
sional setup. To preserve macroscopical averages we keep the centre of mass coordinates
and momentum fixed and work with the six degrees of freedom for the relative momen-
tum and position. Based on the idea of a delocalised particle pair, we assume that the
relative phase space coordinates of the particles are distributed according to a Gaus-
sian. We then choose the new relative position and momentum randomly according
to the Gaussian probability distribution with the additional constraints given by the
symmetries of the problem. Energy and Lz conservation for instance reduce the prob-
lem by two degrees of freedom. If energy and total angular momentum are conserved
the problem is reduced by four degrees of freedom and so on. As the macroscopical
outcome was found to be insensitive towards the details of the microscopical setup we
do not explore this general collisional setup in the present work, but use the simpler
method described above.
8.4 Tests and optimal parameters
How accurately the numerical setup represents the physical picture depends crucially
on the values of the simulation parameters, in particular N˜/N , the time step ∆t and
the smearing widths wr and wp. In all our simulations we use N˜/N = 10, which
is sufficient for temperatures T ≥ 0.2T˜F . The optimal value of ∆t depends on the
physical parameters of the system. Obvious requirements are that the time step must
be smaller than the typical time between two collisions and that the average distance
travelled during a time step must be much smaller than the diameter of the cross section,
〈vrel〉∆t ¿
√
〈σ˜〉/pi. Another constraint is that the time step must be smaller than the
87
half-period with respect to the largest trap frequency, ∆t < pi/ωmax, but for the aspect
ratios considered in this study this condition is much weaker than the other ones. There
is no lower bound on the time step, however the simulation slows down with decreasing
∆t. All tests described below were performed for systems with N = 10000 atoms.
We performed several tests of the simulation to ensure that it accurately captures
the relevant physics. As mentioned above, we checked that total energy conservation
is indeed satisfied for the whole system to a high accuracy. To obtain the correct
dynamical properties, for instance the damping time of excitation modes later on, we
need to ensure that the equilibrium collision rate observed in the simulation corresponds
to the correct theoretical value. The number of collisions in a given time interval can be
very easily obtained from the simulation, simply by counting all test particle collisions
and then multiplying by the ratio (N/N˜). The theoretical value for the collision rate γ
in the presence of Pauli blocking is given by the following integral,
γblock =
∫
d3r
∫ d3ps
(2pi)3
d3ps
(2pi)3
∫
dΩ dσdΩ |ps − ps|fsfs(1− f
′
s)(1− f ′s). (8.34)
After inserting the Fermi-Dirac distribution this integral can be calculated numerically
[87]. Our numerical setup also provides the powerful tool to artificially switch off Pauli
blocking. This allows to separately check for errors related to the general numerical
setup and the calculation of the Pauli blocking factors. Without Pauli blocking the
integral is simpler and can be solved analytically for the Maxwell-Boltzmann distribu-
tion,
γnoblock =
∫
d3r
∫ d3ps
(2pi)3
d3ps
(2pi)3
∫
dΩ dσdΩ |ps − ps|fsfs (8.35)
= N↑N↓ 2ω
3
0
piT 2
(
1 + 1ma2T e
1
ma2T Ei
(
− 1ma2T
))
, (8.36)
where Ei(x) = ∫ x−∞(et/t)dt is the exponential integral. Furthermore we can obtain the
theoretical prediction for the Pauli blocking probability by solving the integral (8.35)
for the Fermi-Dirac distribution. The probability pPauli that two classically colliding
fermions will indeed scatter is then be given by the ratio of γblock to this integral.
To find the optimal value for ∆t for each system we measure the collision rate in the
absence of Pauli blocking for decreasing values of the time step and compare it to the
theoretical prediction. The time step is small enough when we reach good agreement.
It is important to switch off Pauli blocking for the tuning of the time step, as in the
presence of Pauli blocking the collision rate is sensitive to the values of the smearing
widths, which can obscure inaccuracies due to a time step that is too large. After having
found the optimal time step we check that repeated collisions of the same particle pair
are rare. This has always been found to be the case with our collisional setup. We then
88
 0
 1
 2
 3
 4
 5
 0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2
T/T~F
γblock/Nsγnoblock/NspPauli
 0
 2
 4
 6
 8
 10
 0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2
T/T~F
γblock/Nsγnoblock/NspPauli
Figure 8.1: The equilibrium collision rates per particle with and without Pauli blocking, as
well as the Pauli probability for a successful scattering versus temperature for
|k˜Fa| = 1 (left) and |k˜Fa| = 2 (right). The lines correspond to the theoretical
prediction and the symbols to the values obtained with the simulation.
use this optimal value for ∆t to establish the optimal values of the smearing widths. We
first identify the allowed interval for wp and wr given by the conditions (8.31) and (8.32)
and perform the same collision rate matching as described above, this time with Pauli
blocking switched on. We find that the optimal widths lie between wr = wp/m = 1.0lho
for the lowest and wr = wp/m = 2.0lho for the highest temperatures used in our analysis,
where lho = 1/√mω0 is the natural harmonic oscillator length unit.
The measured collision rates for the optimal choice of parameters with and without
Pauli blocking together with the theoretical predictions are shown in Fig. 8.1. For
sufficiently high temperatures the agreement is very good. At very low temperatures it
becomes increasingly difficult to resolve the conditions (8.31) and (8.32) for the Gaussian
smearing widths simultaneously. For larger values of the scattering length this problem
gets worse since larger cross sections can lead to collisions between test particles which
are further apart [87] and therefore the continuous distribution function needs to be
resolved accurately over larger scales.
Another important test is to check that the system thermalises to the correct equi-
librium energy distribution, independently of the initial distribution. In the presence
of Pauli blocking the energy is distributed according to the Fermi-Dirac distribution
(8.8), as shown in section 8.1. Without Pauli blocking the particles will be distributed
according to the Maxwell-Boltzmann distribution (8.20). Figures 8.2 and 8.3 show the
results of this test for a low temperature system with |k˜Fa| = 1 and isotropic trap
frequencies. The parameter values for the time step and the smearing widths are op-
timal. We performed two tests of the thermalisation. First we initialised the system
according to the Fermi-Dirac distribution for T = 0.2T˜F (Fig. 8.2) and ran the sim-
ulation without Pauli blocking. After a short time the system thermalised according
to the Maxwell-Boltzmann distribution for T = 0.31T˜F . Note that the temperatures
of the two distributions are not necessarily equal, since the equipartition theorem does
89
 0
 0.02
 0.04
 0.06
 0.08
 0.1
 0.12
 0  20  40  60  80  100
f s(E
)/N s
E/ω0
start distribution
MB distribution
FD distribution
 0
 0.02
 0.04
 0.06
 0.08
 0.1
 0.12
 0  20  40  60  80  100
f s(E
)/N s
E/ω0
end distribution
MB distribution
FD distribution
 1e-07
 1e-06
 1e-05
 0.0001
 0.001
 0.01
 0  20  40  60  80  100
(ω 03
/E2
)f s(E
)/N s
E/ω0
start distribution
MB distribution
FD distribution
 1e-07
 1e-06
 1e-05
 0.0001
 0.001
 0.01
 0  20  40  60  80  100
(ω 03
/E2
)f s(E
)/N s
E/ω0
end distribution
MB distribution
FD distribution
Figure 8.2: The equilibrium energy distributions without Pauli blocking. The top panel
shows the energy distributions on a linear scale, the bottom panel shows the
energy distributions scaled by ω30/E2 on a logarithmic scale. The start distribu-
tion (left) is Fermi-Dirac and the end distribution (right) is Maxwell-Boltzmann
as expected.
not hold for the Fermi-Dirac distribution. For the Maxwell-Boltzmann distribution we
have
〈E〉 = 2〈Ekin〉 = 3TN, (8.37)
while for the Fermi-Dirac distribution we have
〈E〉 =
∫ ∞
0
(E3/ω30)dE
1 + e(E−µs)/T = −
6T 4
ω30
Li4(−eµs/T ), (8.38)
where Lin(z) =
∑∞
k=1
zk
kn is the polylogarithm. When changing from one distribution to
the other the average energy of the system remains conserved and hence the temperature
of the new equilibrium state is different. Figure 8.3 shows the corresponding results for
the reverse situation: the initial distribution is Maxwell-Boltzmann with T = 0.31T˜F
and Pauli blocking is switched on. It is clearly visible from both figures that the correct
equilibrium distribution is always attained at the end of the simulation. This agreement
improves further at higher temperatures.
Collective excitations emerge when a many-particle system is perturbed away from
equilibrium. Here we will discuss three different excitation modes, the sloshing mode
90
 0
 0.02
 0.04
 0.06
 0.08
 0.1
 0.12
 0  20  40  60  80  100
f s(E
)/N s
E/ω0
start distribution
MB distribution
FD distribution
 0
 0.02
 0.04
 0.06
 0.08
 0.1
 0.12
 0  20  40  60  80  100
f s(E
)/N s
E/ω0
end distribution
MB distribution
FD distribution
 1e-07
 1e-06
 1e-05
 0.0001
 0.001
 0.01
 0  20  40  60  80  100
(ω 03
/E2
)f s(E
)/N s
E/ω0
start distribution
MB distribution
FD distribution
 1e-07
 1e-06
 1e-05
 0.0001
 0.001
 0.01
 0  20  40  60  80  100
(ω 03
/E2
)f s(E
)/N s
E/ω0
end distribution
MB distribution
FD distribution
Figure 8.3: The equilibrium energy distributions with Pauli blocking. The top panel shows
the energy distributions on a linear scale, the bottom panel shows the energy
distributions scaled by ω30/E2 on a logarithmic scale. The start distribution
(left) is Maxwell-Boltzmann and the end distribution (right) is Fermi-Dirac as
expected.
(also known as dipole or Kohn mode), the breathing mode (monopole mode) and the
quadrupole mode. We will confirm that the simulation gives the correct frequencies and
damping properties of these modes. The following tests were performed for a spherical
trap ωx = ωy = ωz = ω0.
The sloshing mode is excited by a small displacement of the centre of mass from
its equilibrium position, or equivalently by a short-lived force represented by an addi-
tional linear term in the potential. It is easy to see that the position and momentum
of the centre of mass of the system must obey the same equations of motion as the
individual particles. As a collision does not change the centre of mass coordinates or
momentum either, 〈r〉 and 〈p〉 must obey the harmonic oscillator equations. Hence
the time evolution of each of the three centre of mass coordinates 〈ri〉 is an undamped
oscillation with the corresponding harmonic oscillator frequency ωi. Figure 8.4 shows
such an oscillation for a system at |k˜Fa| = 1 and T = 0.2T˜F .
The breathing mode can be excited by compressing or expanding the cloud. In a
spherical trap this yields an undamped oscillation of 〈r2〉 with frequency 2ω0 around
its equilibrium value. To see this consider the averages of the kinetic and the potential
91
-0.3
-0.2
-0.1
 0
 0.1
 0.2
 0.3
0 2pi 4pi 6pi 8pi
<
z>
/l z
ωzt
Figure 8.4: Simulation of the equilibrium sloshing mode 〈z〉/lz, where lz = 1/√mωz. In
agreement with the theory the mode is undamped and the oscillation frequency
equals ωz.
energy 〈Ekin〉 = 〈p2〉/2m and 〈Epot〉 = mω20〈r2〉/2. Inserting the harmonic oscillator
equations of motion r˙ = p/m and p˙ = −mω20r we obtain
d
dt(〈Ekin〉 − 〈Epot〉) = −2ω
2
0〈r · p〉, (8.39)
d
dt〈r · p〉 = 2(〈Ekin〉 − 〈Epot〉). (8.40)
These are again the equations for an undamped harmonic oscillator, this time with the
frequency 2ω0. Let us now consider the effect of collisions. Conservation of the total
momentum and the modulus of the relative momentum of the two colliding particles im-
plies that 〈Ekin〉 is unchanged during a collision. By the same reasoning conservation of
the centre of mass coordinates and the relative distance of the two particles implies that
〈Epot〉 is unchanged for an isotropic potential. Since the two particles are approximately
at the same point in space when they collide the term 〈r · p〉 also stays approximately
unchanged due to momentum conservation. Hence the term 〈Ekin〉 − 〈Epot〉 oscillates
around zero with frequency 2ω0 and this oscillation is undamped. This implies that
after a small perturbation 〈r2〉 will perform an undamped oscillation with frequency 2ω0
around its equilibrium value 〈r2〉eq = 〈p2〉/(m2ω20). Figure 8.5 illustrates the breathing
mode for a system with |k˜Fa| = 1 and T = 0.2T˜F .
By perturbing the system via a short-lived small increase in one or several of the
trap frequencies we can excite the quadrupole mode. This is equivalent to a small
momentum shift that equals the gradient of the perturbation potential [87]. We excite
92
 0.98
 0.99
 1
 1.01
 1.02
0 2pi 4pi 6pi 8pi
<
r2
>
/<
r2
>
eq
ω0t
Figure 8.5: Simulation of the normalised equilibrium breathing mode 〈r2〉/〈r2〉eq. In agree-
ment with the theory the mode is undamped and the oscillation frequency equals
2ωz.
the quadrupole mode Q(t) = 〈x2〉 − 〈y2〉 by applying the perturbation ∆px = −cx and
∆py = cy with c = 0.2mω0 in the same way as in [87]. Unlike the sloshing and the
breathing mode this mode is damped. The frequency of the quadrupole mode at high
temperatures approaches the ideal gas value 2ω0, while at low temperatures it is closer
to the hydrodynamic frequency √2ω0 [87]. Figure 8.6 shows the quadrupole mode for
|k˜Fa| = 1 in the high and in the low temperature regime. In both cases the damping
of the mode is clearly visible. The damping is stronger at lower temperatures when
collisions are more frequent. The frequency of the oscillation can be extracted from the
corresponding Fourier transform Q(ω) = ∫∞0 dtQ(t)eiωt and is in agreement with the
theoretical prediction.
93
-0.2
-0.1
 0
 0.1
 0.2
0 2pi 4pi 6pi 8pi
Q(t)
/<r2
> e
q
ω0t
 0
 0.5
 1
 1.5
 2
 0  0.5  1  1.5  2  2.5  3
-
Im(Q
(ω))
ω 0
/<r2
> e
q
ω/ω0
-0.2
-0.1
 0
 0.1
 0.2
0 2pi 4pi 6pi 8pi
Q(t)
/<r2
> e
q
ω0t
 0
 0.5
 1
 1.5
 2
 0  0.5  1  1.5  2  2.5  3
-
Im(Q
(ω))
ω 0
/<r2
> e
q
ω/ω0
Figure 8.6: Left: Simulation of the normalised equilibrium quadrupole mode Q(t)/〈r2〉eq at
high temperature T = 2T˜F (top panel) and at low temperature T = 0.4T˜F (bot-
tom panel). Right: The imaginary part of the corresponding Fourier transforms
of the numerical data gives information about the oscillation frequency.
94
Chapter 9
Collision of two fermionic clouds
We will now apply the numerical setup introduced in the previous chapter to study
the collision of two fermionic clouds with opposite spin polarisation in a cigar shaped
harmonic trap. This research is motivated by the recent experiments of A. Sommer et
al. [37, 40]. Our approach is complementary to the experimental one given our ability
to measure experimentally inaccessible observables such as the local collision rate. A
recent theoretical work using a hydrodynamic approach based on a many-body equation
of state addresses similar issues [90].
The initial distributions are created by sampling a Fermi-Dirac distribution given
by Eq. (8.2). Then each spin distribution is displaced in opposite directions along the
z-axis by d0/2. In the subsequent time evolution, the clouds begin to move towards the
centre under the harmonic trapping force resulting in a collision between them. After
a sufficiently long time the centre of mass energy 2 · 12mω2z(d0/2)2 will be transformed
completely into the internal energy of the gas and a new equilibrium state will be
reached, characterised by a new Fermi-Dirac distribution with temperature Tfinal (and
a chemical potential µfinal). Note that Tfinal is a function of only the atom number,
the initial temperature Tinit and the initial cloud separation d0 and can be calculated
exactly from these values using energy conservation. In the following whenever we refer
to the final temperature of the system we use the theoretical value obtained from the
initial system properties.
Simulations were carried out for a range of initial cloud temperatures 0.2 ≤ Tinit/T˜F ≤
10, interaction strength 0.5 ≤ |k˜Fa| ≤ 10, and initial distances between the centres
of mass of the two clouds 0.4σz ≤ d0 ≤ 16σz, where σz =
√
Tinit/mω2z . All nu-
merical data presented in this chapter was obtained for N = 10000, N˜ = 10N and
ωx/ωz = ωy/ωz = 100.
9.1 Qualitative Behaviour
The behaviour of the clouds during the simulation can be studied by measuring the
distance between their centres of mass d(t) = 〈z↑ − z↓〉(t). As in [37], we find that
before they come to rest at thermal equilibrium, the motion of the clouds exhibits three
typical behaviours, see Fig. 9.1.
95
 0
 1
0 pi 2pi 3pi 4pitωz
d(t)/d0 0.8
 1.2 1.6
 2 γ/ωz
(a) Transmission
 0
 1
0 pi 2pi 3pi 4pitωz
d(t)/d0 0.8
 1.2 1.6
 2 γ/ωz
(b) Intermediate
 0
 1
0 pi 2pi 3pi 4pitωz
d(t)/d0
 20 25
 30 35 γ/ωz
(c) Bounce
Figure 9.1: Bottom panel: The normalised dipole mode d(t)/d0 for the three different be-
haviours: transmission (a), intermediate (b) and bounce (c). Top panel: the
corresponding collision rate per particle γ/ωz measured in the region with
|z| ≤ σz/2 around the trap centre, where σz =
√
Tinit/mω2z .
Transmission: For sufficiently high temperatures and small interactions, the clouds
oscillate through each other (i.e. d(t) crosses zero at short times) with decreasing
amplitude, see Fig. 9.1(a).
Bounce: At low temperatures and strong interactions the clouds bounce off each other
several times (in each bounce the motion of the centre of mass of each cloud is
reversed at short times and without d(t) crossing zero) before a longer period of
slow approach, see Fig. 9.1(c).
Intermediate: Between the transmission and bounce regimes there is a range of tem-
peratures and interactions where the slow approach behaviour is visible from the
start and neither bounces nor transmissions are observed, see Fig. 9.1(b).
The dependence of the behaviour on temperature and interactions is related to the
variation in collision rate γ in the overlap region between the two clouds. As the collision
rate decreases, the system behaviour changes from the bounce regime, to intermediate
and finally to the transmission regime. From the top panel of Fig. 9.1 we see that,
in the bounce regime, the oscillations in the collision rate integrated over a volume
in the overlap region follow closely the oscillations of d(t), whereas no such variation
is apparent in the transmission regime. In addition, we can compare the collision
rate with the typical timescale for macroscopic motion (ω−1z ). We see that the gas is
strongly hydrodynamic in the bounce regime (ωz/γ ¿ 1) and becoming collisionless in
the transmission regime (ωz/γ ∼ 1).
Note that no artificial repulsion between the atoms was included in the numerical
setup. The repulsive behaviour of the centres of mass of the clouds in the bounce
regime arises naturally from the purely attractive interaction between the individual
96
 0
 20
 40
 60
 80
 0  0.5  1  1.5  2  2.5  3  3.5  4  4.5
(T/
T F
) fin
al
|kFa|final
Figure 9.2: The transition between bounce and intermediate regimes (filled symbols, lines
to guide the eye) and between intermediate and transmission regimes (empty
symbols). Red circles correspond to d0 = 43.1lz, blue triangles to d0 = 64.6lz
and green squares to d0 = 129.3lz, where lz = 1/√mωz. It is clearly visible that
the intermediate-transmission transition is independent of d0. The dashed line
corresponds to constant relaxation time 1/τdip = 1.83ωz.
particles. We have also seen that the bounce can be understood purely in terms of
semiclassical collisions, without the need for mean fields or other more complicated
many-body effects. In particular our approach predicts that all quantities considered
here depend only on the square of the scattering length and not on its sign, a fact
that was also observed in [37]. In contrast, in [90] the bounce is understood from the
equation of state of the so-called “repulsive” branch and therefore, according to that
work, would in principle exhibit different behaviour on opposite sides of the resonance.
9.2 Transitions between regimes
Figure 9.2 shows the transitions between the three regimes of behaviour in the (T/TF )final,
|kFa|final plane for different values of d0 in units lz = 1/√mωz. Remember that as the
density is a function of the temperature, TF and kF defined in (8.12) change during the
simulation and hence kFa varies during the evolution. The quantities given here are
equilibrium values for t →∞.
The transition between transmission and intermediate regimes is defined as d(t)
reaching but not crossing zero at short times. It is clearly visible from Fig. 9.2 that
this transition is independent of the initial separation d0 and hence can be understood
entirely from the final equilibrium properties of the system. More precisely, it can be
97
understood as a consequence of the change in the relaxation time τdip of the spin dipole
mode of the d0 = 0 system, which is closely related to the collision rate per atom. For
sufficiently high temperatures T & T˜F the relaxation time can be calculated from the
Maxwell-Boltzmann distribution and equals
1
τdip =
2N
3piT 2ω
3
0f
( 1
ma2T
)
, (9.1)
where f(y) = 1 − y + y2eyΓ(0, y), as was shown in [91]. The incomplete Γ-function is
hereby defined as Γ(0, y) = −Ei(−z) = ∫∞x (e−t/t)dt. By comparison with (8.36) we
can see that the relaxation time is closely related to the equilibrium collision rate for
the Maxwell-Boltzmann distribution. For sufficiently low collision rate (ωzτdip À 1),
the gas can be said to be collisionless and therefore the clouds undergo independent
oscillations without interacting strongly with each other. In the transmission regime
d(t) is simply the solution of the classical damped harmonic oscillator equations [91],
d(t) = e−t/2τdip(A sin(ωt) +B cos(ωt)), (9.2)
where A and B are given by the initial conditions and the oscillation frequency equals
ω = √ω2z − 1/(2τdip)2. Hence the oscillation frequency is related to the damping. When
1/τdip → 0 the frequency becomes exactly ωz and there is no damping (non-interacting
gas). At finite τdip the frequency is smaller than ωz. For 1/τdip > 2ωz the dipole mode
is overdamped [91, 92] as ω becomes ill-defined. To compare with the simulations, we
can calculate the value of τdip for the various points lying on the curve separating the
transmission and intermediate regimes of Fig. 9.2 using Eq. (9.1). We find that they
lie on the curve of constant 1/τdip = 1.83(1)ωz.
The transition between the intermediate and the bounce regime is defined to occur
when the first bounce ceases to reverse the motion of the clouds, or in other words
when d(t) ceases to have a minimum and becomes a monotonically decreasing function
of t. This transition depends on d0. In the bounce regime we typically see an initial
strong collision followed by oscillations of d(t) which eventually die out as d(t) → 0.
This oscillatory behaviour continues into the intermediate regime. As the damped
harmonic oscillator can never produce such a bounce solution, this regime must be the
consequence of a different mechanism. In the following section we will relate it to a
non-linear coupling between the spin dipole mode and the breathing mode.
9.3 Coupling between excitation modes
We have already discussed collective excitations of systems close to equilibrium in sec-
tion 8.4. The spin dipole mode is an oscillation of the centre of mass of the system. In
our case we consider d(t)/d0, the distance between the z-coordinates of the centres of
98
 0
 1
0 pi 2pi 3pi 4pitωz
d(t)/d0<z2>/<z2>eq
(a) Transmission
 0
 1
0 pi 2pi 3pi 4pitωz
d(t)/d0<z2>/<z2>eq
(b) Intermediate
 0
 1
 2
0 pi 2pi 3pi 4pitωz
d(t)/d0<z2>/<z2>eq
(c) Bounce
Figure 9.3: The normalised dipole mode d(t)/d0 (red solid lines) and breathing mode
b(t)/b∞ (blue dashed lines) for the three different behaviours: transmission (a),
intermediate (b) and bounce (c).
mass of the two clouds, normalised by their initial separation. We will compare it to
the axial breathing mode b(t) = 〈z2↑+z2↓〉(t). Figure 9.3 shows plots in the three regimes
of the normalised amplitude of the dipole mode d(t)/d0 and normalised amplitude of
the breathing mode b(t)/b∞. We can see from Fig. 9.3(c) that in the bounce regime
the frequency of d(t) is identical to the frequency of b(t). This suggests the existence
of a non-linear coupling between the two modes. As we move towards the intermediate
regime shown in Fig. 9.3(b) the frequency of the dipole mode becomes ill-defined due
to the strong damping. In the transmission regime, see Fig. 9.3(a), it becomes closer
to the frequency of the dipole mode of the non-interacting gas ωz, as discussed in the
previous section. This indicates that the dipole mode decouples from the breathing
mode.
For a more quantitative analysis we fit d(t) with the function
d(t) = Be−t/C (1 +De−t/E sin (ωt+ φ)) . (9.3)
The first term is related to the spin drag coefficient measured in [37], which we will
analyse in section 9.4. It dominates the overdamped behaviour of d(t) at long times
with a characteristic timescale C. The second term was not quantitatively analysed
in [37] since it takes into account the short time behaviour which includes a damped
oscillation with a typically shorter damping time scale E.
For displacements d0 & 2σz the values of the fit parameters depend significantly
on the range of the data included in the fit. To illustrate this Fig. 9.4 shows fits for
different data ranges for the system with Tinit = 0.4T˜F , |k˜Fa| = 1 and d0 = 3σz.
Table 9.1 compares the corresponding fit parameter values. The frequency ω/ωz for
instance becomes smaller with increasing cut-off times. Since the frequency shows this
weak time dependence we ignore early times by imposing a cut-off dc on the amplitude
and fitting only times t > tc for which d(t) < dc. Since the function d(t) is not monotonic
99
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1
 1.1
0 pi 2pi 3pi 4pi 5pi 6pi
d(t)/
d 0
tωz
data
fit
(a) full data range, tcωz = 0
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1
 1.1
0 pi 2pi 3pi 4pi 5pi 6pi
d(t)/
d 0
tωz
data
fit
(b) without first bounce, tcωz ≈ 1.1pi
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1
 1.1
0 pi 2pi 3pi 4pi 5pi 6pi
d(t)/
d 0
tωz
data
fit
(c) without first two bounces, tcωz ≈ 2.3pi
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1
 1.1
0 pi 2pi 3pi 4pi 5pi 6pi
d(t)/
d 0
tωz
data
fit
(d) without first three bounces, tcωz ≈ 3.5pi
Figure 9.4: Dependence of the best fit of the normalised dipole mode d(t)/d0 to Eq. (9.3)
on the range of the data used for the fit. All data were obtained for the system
Tinit = 0.4T˜F , |k˜Fa| = 1, d0 = 3σz.
tcωz B/d0 Cωz D Eωz ω/ωz φ
full data range 0 0.668 32.7 0.488 6.59 1.71 0.459pi
without first bounce 1.1pi 0.652 34.9 0.425 8.11 1.64 0.606pi
without first two bounces 2.3pi 0.652 34.9 0.392 8.66 1.62 0.696pi
without first three bounces 3.5pi 0.655 34.5 0.379 8.86 1.59 0.806pi
Table 9.1: Dependence of the best fit parameters to Eq. (9.3) on the data range used for the
fit. All data were obtained for the system Tinit = 0.4T˜F , |k˜Fa| = 1, d0 = 3σz.
the corresponding time cut-off tc is not continuous with changing d0, which leads to a
small systematic error. We estimate this error from the cut-off dependence of the fit
parameters. The statistical error of the fit is several orders of magnitude smaller and
hence negligible.
100
 0
 0.5
 1
 1.5
 2
 2.5
 3
 0  2  4  6  8  10
D’/2
D
d0/σz
 0.4
 0.6
 0.8
 1
 1.2
 1.4
 1.6
 0  2  4  6  8  10
E’/E
d0/σz
 0.9
 0.95
 1
 1.05
 1.1
 0  2  4  6  8  10
ω’
/ω
d0/σz
 0
 0.5
 1
 1.5
 2
 0  2  4  6  8  10
φ’/φ
d0/σz
Figure 9.5: Comparison between the fit parameters in (9.3) and (9.4) for |k˜Fa| = 1, Tinit =
0.4T˜F and different values of initial separation. The numerical results confirm
the relations (9.5).
The oscillations of the axial size of the cloud are fitted using the dependence
b(t)/b∞ = 1 +D′e−t/E′ sin (ω′t+ φ′) . (9.4)
From the fit to our simulations, we observe that
D′ ' 2D, (9.5a)
E ′ ' E, (9.5b)
ω′ ' ω, (9.5c)
φ′ ' φ, (9.5d)
as shown in Fig. 9.5. This result can be explained by noting that in the hydrodynamic
and the collisionless regimes, the spin dipole and breathing modes can be described by
a shift and a scaling of the phase space density. We assume that this also applies here
and take for the variation of the density of the spin species s along the z-axis the ansatz
ns(z, t) = αns(αz ± β, 0), (9.6)
101
where α = α(t) and β = β(t) represent the breathing and spin dipole modes respectively.
Using this ansatz, we can express b(t) and d(t) as functions of α and β:
〈z〉 =
∫
zns(z, t)dz, (9.7)
=
∫
zαns(αz ± β, 0)dz, (9.8)
= 1α
∫
(u± ∓ β)ns(u±, 0)du±, with u± = αz ± β, (9.9)
= ∓βα, (9.10)
using ∫ ns(u)du = 1 and
∫ uns(u)du = 0. Hence
d(t) = 〈z↑ − z↓〉(t) = 2β(t)α(t) . (9.11)
Similarly
〈z2〉 =
∫
z2αns(αz ± β, 0)dz, (9.12)
= 1α2
∫
(u± ∓ β)2ns(u±, 0)du±, (9.13)
= 1α2
(〈u2〉eq. + β2
) , (9.14)
as the linear term again vanishes. If we also ignore the second order contribution in β
we obtain
b(t)
b∞ =
〈z2↑ + z2↓〉(t)
〈z2〉eq. ≈
1
α2(t) . (9.15)
Now we assume that β(t) is overdamped,
β(t) = Be−t/C , (9.16)
and that α(t) is a damped oscillation,
α(t) = 1 +De−t/E sin(ωt+ φ), (9.17)
and then by expanding Eqs. (9.11) and (9.15) to leading order we obtain the time
dependence in Eqs. (9.3) and (9.4) together with the corresponding relations between
the coefficients (9.5).
We are also able to study the temperature dependence of the frequency of the two
modes, see Fig. 9.6. We obtain the frequency by fitting to the functions (9.3) for the
dipole mode and (9.4) for the breathing mode. As we increase the temperature, the
frequency of the spin dipole and breathing modes increases from
√
12/5ωz ≈ 1.55ωz,
102
 1.5
 1.55
 1.6
 1.65
 1.7
 1.75
 1.8
 1.85
 1.9
 1.95
 2
 2.05
 0  5  10  15  20  25
ω
/ω
z
(T/TF)final
Figure 9.6: The frequency ω/ωz of the dipole mode d(t) (red circles) and the breathing mode
b(t) (blue triangles) versus the final temperature for |k˜Fa| = 1. All data were
obtained for equal initial temperature Tinit = 0.4T˜F by varying d0. The solid
line is the prediction from [93].
the frequency of the hydrodynamic axial breathing mode [93, 94], and approaches 2ωz,
the non-interacting value. This is in agreement with the value for the bounce frequency
at unitarity ω = 1.63(2)ωz measured in [37]. At higher temperatures the spin dipole
mode frequency becomes ill-defined due to large damping. If we continue to increase the
temperature, the damping becomes progressively smaller and the dipole mode frequency
approaches ωz, the value for the non-interacting gas, as was discussed above. We
also compare our results with an earlier prediction for the frequency of the breathing
mode [93]. The two dependencies are close to each other although some discrepancy
remains. We attribute it to the fact that to estimate the collisional rate in the cloud
Ref. [93] neglects the Pauli principle and assumes a Gaussian phase-space density. The
temperature dependence of the frequency remains the same if we change the scattering
length or the initial temperature of the system, as shown in Fig. 9.7.
9.4 Transport coefficients
Transport coefficients characterise the macroscopic response of a system towards an
external perturbation. Usually we consider linear response to weak perturbations with
the transport coefficient being the proportionality factor. In this section we will discuss
how different transport coefficients can be calculated from the numerical data.
The spin drag coefficient Γsd characterises the linear response in the relative spin
103
 1.5
 1.6
 1.7
 1.8
 1.9
 2
 2.1
 0  5  10  15  20  25
ω/ω
z
(T/TF)final
|k~Fa| = 1|k~Fa| = 2
 1.4
 1.5
 1.6
 1.7
 1.8
 1.9
 2
 2.1
 2.2
 2.3
 5  10  15  20  25
ω/ω
z
(T/TF)final
Tinit = 0.4T~FTinit = 1.25T~F
Figure 9.7: The frequency ω/ωz of the dipole mode d(t) versus the final temperature for
Tinit = 0.4T˜F and two different values of the scattering length (left), and for
|k˜Fa| = 1 and two different values of the initial temperature (right). No depen-
dence on either the scattering length or the initial temperature is visible.
current towards a weak external force. It is defined through
〈F↑↓〉 = −Γsd〈p↑ − p↓〉, (9.18)
where 〈F↑↓〉 = −(∇V (〈r↑〉) − ∇V (〈r↓〉)) is the force between the centres of mass of
the two fermionic clouds [95]. If the clouds are displaced along the z-direction in a
harmonic external potential we obtain
ω2zd+ Γsdd˙ = 0, (9.19)
for small displacements from equilibrium. The solution to this equation is
d(t) ∝ e−tω2z/Γsd , (9.20)
in agreement with the overdamped behaviour (9.16) near equilibrium which is seen in
the experiment and the numerical simulation at sufficiently long times. Hence the spin
drag coefficient can be obtained from the fit coefficient C via
Γsd = Cω2z . (9.21)
An explicit form for the spin drag of the homogeneous gas at non-degenerate tempera-
tures can be obtained by solving the Boltzmann equation in the classical regime with a
variational approach, as shown for instance in the supplementary material to [37] and
references therein. One obtains
Γ(hom)sd =
8
3nσ(a, T )
( T
pim
)1/2
, (9.22)
104
where
σ(a, T ) = 4pia2
∫
du u
5e−u2
1 + u2ma2T (9.23)
= 2pimT
(
1− 1ma2T −
1
(ma2T )2 e
1
ma2T Ei
(
− 1ma2T
))
. (9.24)
For weak interactions the cross section stays nearly constant, σ(ma2T ¿ 1) = 4pia2,
and the spin drag is proportional to the square root of the temperature. In the limit of
large scattering length (ma2T À 1) expression (9.24) simplifies to
σ(ma2T À 1) = 2pimT , (9.25)
which is the square of the thermal de Broglie wavelength. Therefore the dimensionless
spin drag coefficient at unitarity is
Γ(hom)sd
εF =
32√2
9pi3/2
( T
TF
)−1/2
≈ 0.90
( T
TF
)−1/2
. (9.26)
At very low temperatures there is no reliable microscopic theory that could predict the
spin drag coefficient. In [38] the spin relaxation time, which is inversely proportional to
the spin drag, is calculated also at low temperature using Fermi liquid theory to include
the effect of strong correlations. The calculation reveals that Γ(hom)sd has a maximum in
the strong coupling limit. At unitarity this maximum is found at temperatures around
the Fermi temperature TF . At low temperatures the spin drag behaves according to
Γ(hom)sd /εF ∼ (T/TF )2 at unitarity, and Γ(hom)sd /εF ∼ (T/TF )2(kFa)2 for weak interac-
tions. This is in agreement with the experimental observations [37]. To understand the
temperature dependence of the spin drag qualitatively note that for a classical gas the
cross section increases with decreasing temperature. The maximum of the spin drag
can be understood as a consequence of the onset of Pauli blocking at low temperatures,
which causes the cross section to decrease again.
In the presence of a non-uniform external potential the spin drag will depend on
the spatial position. From a fit to d(t) we can obtain only the trap average Γsd which
is related to the local spin drag Γsd(r) via
Γsdd˙ = Γsd〈vrel,z(r)〉 = 〈Γsd(r)vrel,z(r)〉. (9.27)
The z-component of the relative velocity vrel,z consists of the equilibrium term and
the drift velocity term ±vd(r) describing the directed motion of the spin up and spin
down atoms towards the trap centre in response to the displacement of the clouds. The
average equilibrium velocity is zero, due to symmetry of the equilibrium distribution.
105
 0
 0.02
 0.04
 0.06
 0.08
 0.1
 0  5  10  15  20  25
Γ sd
/ε F
(T/TF)final
Tinit = 0.4T~FTinit = 1.25T~F
 0
 0.05
 0.1
 0.15
 0.2
 0  5  10  15  20  25
Γ sd
/ε F
(T/TF)final
Tinit = 0.4T~FTinit = 1.25T~F
Figure 9.8: The average spin drag Γsd versus final temperature for |k˜Fa| = 1 (left) and
|k˜Fa| = 2 (right). Red circles denote simulation results obtained for Tinit =
0.4T˜F and blue triangles simulation results obtained for Tinit = 1.25T˜F , by
varying the initial separation. The lines are fits of the data for Tinit = 1.25T˜F
to Γ(hom)sd /α, where the homogeneous spin drag is given by (9.22) and α is a free
fit parameter.
Hence we can write for the average spin drag,
Γsd = 〈Γsd(r)vrel,z(r)〉〈vrel,z(r)〉 =
∫ d3rΓsd(r)n(r)vd(r)∫ d3rn(r)vd(r) . (9.28)
The local spin drag Γsd(r) can be approximated by the homogeneous value (9.22) with
the uniform density n replaced by the local density n(r). This yields
Γsd = Γsd(0)
∫ d3rn2(r)vd(r)
n(0) ∫ d3rn(r)vd(r) . (9.29)
Hence the average spin drag coefficient is proportional to the local value at the trap
centre, Γsd = Γsd(0)/α, where the coefficient α is given by
α = n(0)
∫ d3rn(r)vd(r)∫ d3rn2(r)vd(r) (9.30)
and can only be calculated by making an assumption about the drift velocity profile of
the gas [37]. In a harmonic potential, one obtains α = 23/2 ≈ 2.83 assuming a uniform
drift velocity, and α = 25/2 ≈ 5.66 if one assumes a quadratic drift velocity profile.
Compared to the local spin drag at the trap centre the average spin drag is lower by
the factor α.
Our results for the spin drag coefficient for |k˜Fa| = 1 and |k˜Fa| = 2 are presented
in Fig. 9.8 together with a fit to the analytical prediction (9.22) at high temperatures.
A key observation is that, unlike the oscillation frequency, the average spin drag shows
a dependence on the initial temperature and displacement and is not purely a function
of the final temperature. We compare two data sets for each scattering length, one
106
obtained at a low initial temperature Tinit = 0.4T˜F and one at a high initial temperature
Tinit = 1.25T˜F . The data at the high initial temperature is fitted well by the analytical
prediction for Γ(hom)sd /α, where the homogeneous spin drag is given by (9.22) and α is
a free fit parameter. For |k˜Fa| = 1 the fit yields α = 3.34(8) and for |k˜Fa| = 2 the fit
yields α = 3.52(11). These two values are consistent with each other within error. If
we fit both data sets together we obtain α = 3.35(15).
The high temperature tail of the data obtained using Tinit = 0.4T˜F and large initial
separations is on the other hand not well-fitted by the theoretical prediction. We
attribute this to the fact that for large initial separations non-linear effects strongly
influence the drift velocity profile and hence the points in the tail do not all correspond
to the same value of α. This obscures the proportionality relation between the average
spin drag and the homogeneous value. In support of this argument note that the points
for Tinit = 0.4T˜F and final temperatures between approximately TF and 5TF lie close to
the fit to the analytical prediction. Only when the final temperature increases further,
for which larger initial separations are needed, does the data start to deviate.
We can also see the predicted maximum of the spin drag around the Fermi temper-
ature. The spin drag then starts to decrease again as we go to even lower temperatures
but this regime is not yet fully accessible with the simulation.
The values for α extracted from the fits lie between the predictions for a uniform and
a quadratic velocity profile. In [37] a higher value for α was found at unitarity, closer
to the prediction for a quadratic drift velocity profile. In the experimental setup fast
particles evaporate from the trap, such that the time-evolution of the clouds is isother-
mal. In our numerical setup the particle number is conserved, but not the temperature.
This affects the drift velocity profile and hence the relation between the average and
the local spin drag coefficients. In our setup the average spin drag coefficient is closer
to the homogeneous value.
It is possible to access also other transport coefficients, for instance the spin diffusion
coefficient Ds which is given by Fick’s law,
〈j↑↓〉 = −Ds∇(n↑ − n↓), (9.31)
where 〈j↑↓〉 = n〈v↑ − v↓〉 is the spin current [37, 38]. Note that in the presence of an
external potential this definition is in terms of the local densities. Using values at the
trap centre the spin diffusion coefficient can be extracted from the simulation data via
Ds = −n(0)d˙g(0) =
n(0)ω2zd
Γsdg(0) , (9.32)
where g(z) = ∂z(n↑ − n↓) is the slope of the local density difference at x = y = 0.
The above formula is only valid close to equilibrium when d(t) decreases exponentially
with time satisfying Eq. (9.19). The quantities d, n(0) and g(0) are time-dependent,
107
however these dependencies should cancel, such that Ds is constant close to equilibrium.
Similarly one can also obtain the spin susceptibility via
χs = g(0)mdω2z
. (9.33)
The calculation of these and other transport coefficients will be subject of future re-
search.
108
Chapter 10
Conclusions and Outlook
In the first part of this thesis we discussed a Monte Carlo calculation of the critical
temperature and other thermodynamic observables of the homogeneous unitary Fermi
gas with equal and unequal chemical potentials for the two spin components. This work
contains two main improvements compared to previous DDMC studies: a new update
scheme which is less susceptible to autocorrelations than the original worm algorithm,
and an alternative procedure for extracting the thermodynamic limit of the critical
temperature from the numerical data that avoids a systematic error present in several
other studies. Due to these improvements we could provide a more accurate estimate
on the critical temperature than the previously established value. We also generalised
the DDMC algorithm to the imbalanced case using the sign quenched method allow-
ing for the first time to numerically access the thermodynamic properties of the spin
imbalanced unitary Fermi gas.
Our most precise values (using data at both zero and non-zero imbalance) for the
temperature, the chemical potential, the energy per particle and the contact density of
the balanced unitary Fermi gas at the critical point are
Tc/εF = 0.171(5), (10.1)
µ/εF = 0.429(7), (10.2)
E/EFG = 0.440(15), (10.3)
C/ε2F = 0.1101(9). (10.4)
The value for the contact density was the first numerically obtained finite temperature
result for this quantity.
For unequal chemical potentials in the two components we extracted the dependence
of the critical temperature on the imbalance h = ∆µ/εF close to the balanced limit.
Our analysis is consistent with Tc/εF = const. for the range of imbalances considered
(h . 0.2). We also derived a lower bound on the leading order term in the expansion
of the critical temperature Tc(h)− Tc(0) > −0.5εFh2. With the additional assumption
that the linear dependence of the lattice critical temperature Tc/εF on ν1/3 remains
unchanged in the presence of a small imbalance, a tighter lower bound of Tc(h) −
Tc(0) > −0.04εFh2 could be obtained. We also analysed the behaviour of the energy,
109
the chemical potential and the contact density in the presence of an imbalance. In all
cases we saw good agreement with experimental results.
We also studied the above mentioned quantities for the balanced gas at temperatures
above and below Tc. We found that the chemical potential decreases with increasing
temperature, while the energy and the contact density increase with increasing temper-
ature. This is still work in progress; in particular we expect to increase the precision
of our results when data at the lowest filling factor used for the analysis at Tc becomes
available. So far we have not studied the spin imbalanced gas at temperatures beyond
Tc, which will be explored in future work.
Our DDMC study could provide results for several thermodynamic observables with
an accuracy of a few percent. Nevertheless there is still much room for improvement.
A major source of statistical errors is the continuum extrapolation of the data. Due
to lattice discretisation effects observables calculated at finite values of the chemical
potential do not equal the physical observables in the continuum limit. Since due
to technical limitations we cannot use arbitrarily large lattice sizes, we cannot work
arbitrarily close to continuum. The slope of the continuum extrapolation is large due
to the significant difference between the lattice dispersion relation and the parabolic
continuous dispersion relation. Reducing this slope would allow us to get more precise
continuum values without having to go to very low values of the chemical potential.
To achieve this a more complex dispersion relation than the one used in this work
can be invoked. We have already briefly discussed dispersion relations that completely
suppress leading order lattice corrections and aim to investigate such dispersion relations
in future work.
The efficiency of the simulation can be further increased by inventing more elaborate
update schemes. For instance updates that insert or remove entire vertex chains into
the configuration would change the configuration significantly in only one Monte Carlo
step. This would further reduce errors associated with autocorrelations and speed the
simulations.
The improved DDMC algorithm with sign quenching also offers the intriguing pos-
sibility to explore the case of unequal masses of the two species. If m↑ 6= m↓ the
dispersion relations are different for the two components and the mass ratio enters as a
new parameter. As in the spin-imbalanced case the two matrix determinants no longer
need to be identical, so that sign quenching is required. The mass ratio is expected to
influence the behaviour of the system significantly, so that exploring the phase diagram
promises many new interesting insights.
To access dynamical properties of a Fermi gas above the critical temperature and at
finite values of the scattering length, we have implemented and tested a numerical sim-
ulation of the Boltzmann equation. With this simulation we have studied the collision
of two spin polarised Fermi clouds in a cigar shaped harmonic trap. In agreement with
recent experimental observations [37] we identified three distinct regimes of behaviour
110
which can be characterised in terms of the centre of mass separation between the two
clouds. For weak interactions the clouds pass through each other. If interactions are in-
creased they approach each other exponentially and for strong interactions they bounce
off each other several times. We characterised the transitions between these regimes as
a function of interaction strength and temperature, and related them to the collision
rate in the overlap region between the clouds. We then explained the occurrence of the
bouncing as the consequence of a non-linear coupling between the spin dipole mode and
the axial breathing mode, which is enforced by collisions. Such a non-linear coupling
between collective modes has, to our knowledge, never before been studied in a Fermi
gas. We also determined the frequency of the bounce as a function of the final temper-
ature of the equilibrated system. Finally, we discussed how spin transport coefficients
can be extracted from the numerical data. We calculated the trap average of the spin
drag coefficient and found good agreement with the analytical prediction.
In future we aim to extend this study to other closely related problems such as
the collision of clouds with unequal populations [40], between clouds of atoms with
unequal masses [96], and to the Fermi liquid regime at T ¿ TF , using the Landau-
Boltzmann equation. The low temperature regime is especially fascinating as, although
it is experimentally accessible, very few analytical results are known.
Our numerical setup gives us precise control over all system parameters and allows
us to monitor quantities not easily accessible in an experiment. We can follow the
trajectory of any given particle and obtain local collision rates, density and velocity
profiles as a function of time. In particular the velocity profile is interesting as it
gives insight about the relation between homogeneous and trap averaged transport
coefficients. Without knowledge of the drift velocity no theory for the average spin
drag in a trap can be formulated. We plan to resolve this question by a detailed
study of the time-dependent drift velocity profile as a function of the temperature, the
scattering length and the initial separation of the clouds. Apart from the spin drag we
intend to calculate also other transport coefficients like the spin diffusion coefficient or
the spin susceptibility.
So far we have not studied the dynamics of the trapped unitary Fermi gas. At
unitarity the s-wave cross section diverges for small relative momenta of the two collid-
ing atoms. This regime regime can be made accessible by enforcing a low-momentum
cut-off, which we aim to explore this in future work.
111
Appendix A
Matrix update procedure
In [50] Rubtsov et al. develop a fast update algorithm designed for the types of matrix
updates described in chapter 5. It is based on calculating the inverse of the matrices,
which gives the determinants as a side product. The calculation of the determinant of
an arbitrary N×N matrix requires order N3 operations. If a matrix, which determinant
is already known, is modified, e.g. by changing one or several rows and/or columns, the
calculation of the determinant of the updated matrix can sometimes be performed in
less than N3 operations. In the following, we will derive all three kinds of update and
give the number of operations needed to perform them. Note that Einstein’s summation
convention is not used: all sums are indicated explicitly and a repeated index is not
necessarily summed over.
A.1 The Sherman-Morrison Formula
Suppose that we have already inverted a matrix A and want to obtain the inverse of a
modified matrix A′, that differs from A in some elements, e.g. in a row or a column. A
matrix inversion from scratch requires order N3 operations, but having already inverted
A we can hope to calculate the inverse of A′ more efficiently. If the change in the matrix
is of the form
A′ = A+ u⊗ v, (A.1)
where u⊗ v is a matrix whose (i, j) element is the product uivj, the Sherman-Morrison
formula gives the inverse of A′ as
(A+ u⊗ v)−1 = A−1 − (A
−1u)⊗ (v · A−1)
1 + v · A−1u . (A.2)
112
It can be derived as follows [97]:
(A+ u⊗ v)−1 = (1 + A−1 · u⊗ v)−1 · A−1 (A.3)
= (1− A−1 · u⊗ v + A−1 · u⊗ v · A−1 · u⊗ v ∓ . . .) · A−1 (A.4)
= A−1 − A−1 · u⊗ v · A−1 + A−1 · u⊗ v · A−1 · u⊗ v · A−1 ∓ . . .(A.5)
= A−1 − A−1 · u⊗ v · A−1(1− Λ + Λ2 ∓ . . .) (A.6)
= A−1 − (A
−1 · u)⊗ (v · A−1)
1 + Λ , (A.7)
where Λ ≡ v · (A−1u) and Taylor expansion was used in (A.4). In (A.6) the scalars Λ
were factored out and in (A.7) the series was written as (1+Λ)−1 using Taylor expansion
again.
Let us count the number of operations required for the calculation of (A+u⊗ v)−1,
first for arbitrary vectors u and v and then for the special cases u = en or v = en,
which correspond to changing all elements of the nth row or nth column of the matrix
respectively. The quantities that need to be computed are
A−1u (vector) → N2 mult., (N − 1)2 add.
v · A−1 (vector) → N2 mult., (N − 1)2 add.
Λ = (v · A−1)u (scalar) → N mult., N − 1 add.
A−1u
1+Λ (vector) → N mult., 1 add. (1 + Λ)
A−1u
1+Λ ⊗ (v · A−1) (matrix) → N2 mult., (1 for each matrix element)
A−1 − A−1u1+Λ ⊗ (v · A−1) (matrix) → N2 add., (1 for each matrix element)
So, in general, this method requires order 3N2 multiplications and 3N2 additions. Mod-
ifying a single row or column of A corresponds to setting either u or v to a unit vector
en. In this case the computation of either A−1u or v · A−1 reduces to just picking the
nth column or the nth row of A−1 respectively. This saves us N2 multiplications and
(N −1)2 additions, so that the total cost reduces to order 2N2 multiplications and 2N2
additions. Unfortunately, a change in one row and one column of a matrix cannot be
written in the form u⊗ v, so that for this update we need to perform the row and the
column update separately, yielding order 4N2 multiplications and additions.
The worm type updates involve additions and removals of one row and one column.
This can be reformulated in terms of changing an existing row and column in the
following way: We can first add a “unit” row n and a “unit” column m to A, with
Anj = δmj and Aim = δin. This does not change the determinant. The inverse matrix is
also almost unchanged, it only gets an extra “unit” row m and “unit” column n, with
A−1mj = δnj and A−1in = δim (note that m and n get swapped for the inverse). Now we
can modify the added row and column in the usual way using the Sherman-Morrison
formula. Removing a row and a column can be done analogously: We change the row
113
n and column m we want to remove to Anj = δmj and Aim = δin using Sherman-
Morrison and then drop them and also drop the corresponding row m and column n
in the inverse matrix. The number of operations for this type of updates is also order
4N2 multiplications and same for additions.
The updates proposed in [50] are an extension of the Sherman-Morrison formula.
In the following three sections, we will derive the three kinds of update and calculate
their complexity.
A.2 Adding one row and one column
Let us first consider the update that adds an (N + 1)st row and column to an N ×N
matrix A (the order of rows and columns is not important, since a permutation of
rows/columns of A yields a permutation of columns/rows of A−1). As described at the
end of the previous section, we first extend the N×N matrix A to an (N +1)× (N +1)
matrix with AN+1,j = Ai,N+1 = 0 for i, j = 1, . . . , N and AN+1,N+1 = 1,
A =


A1,1 · · · A1,N
... ...
AN,1 · · · AN,N

→


0
A ...
0
0 · · · 0 1


=


A1,1 · · · A1,N 0
... ... ...
AN,1 · · · AN,N 0
0 · · · 0 1


.
(A.8)
This does not change the determinant and the inverse of A also gets extended in the
same way,
A−1 =


A−11,1 · · · A−11,N... ...
A−1N,1 · · · A−1N,N

→


0
A−1 ...
0
0 · · · 0 1


=


A−11,1 · · · A−11,N 0... ... ...
A−1N,1 · · · A−1N,N 0
0 · · · 0 1


.
(A.9)
The extended matrices will be denoted with the same letters as the original ones, since
from here on only the former will be used. The extension enables us to perform matrix
additions and multiplications of the original and the updated matrix, for it is only
possible to perform these operations with matrices of the same size. We want to find
the inverse of the updated matrix A′ = A+∆, where
∆ =


0 · · · 0 A1,N+1
... ... ...
0 · · · 0 AN,N+1
AN+1,1 · · · AN+1,N AN+1,N+1 − 1


. (A.10)
114
As in (A.3) we write
A′−1 = (A+∆)−1 = A−1(1 + ∆A−1)−1. (A.11)
The next step in the derivation of the Sherman-Morrison formula was Taylor expanding
the inverse of the term in brackets. In our case however, since ∆ cannot be decomposed
into a direct product of two vectors, this will not yield the desired result. Instead we
can try to invert 1 + ∆A−1 directly. For this consider
∆A−1 =


0 · · · 0 A1,N+1
... ... ...
0 · · · 0 AN,N+1
AN+1,1 · · · AN+1,N AN+1,N+1 − 1


·


A−11,1 · · · A−11,N 0... ... ...
A−1N,1 · · · A−1N,N 0
0 · · · 0 1


=


0 · · · 0 A1,N+1
... ... ...
0 · · · 0 AN,N+1
R1 · · · RN AN+1,N+1 − 1


, (A.12)
where Rj ≡
∑
k AN+1,kA−1k,j. Thus
1 + ∆A−1 =


1 0 · · · 0 A1,N+1
0 1 0
... . . . ... ...
0 · · · 1 AN,N+1
R1 · · · RN AN+1,N+1


. (A.13)
This matrix is naturally partitioned into four blocks, an N × N block in the upper
left corner, containing the unit matrix, a 1 × 1 “block” containing the single element
AN+1,N+1 and the two remaining blocks (1×N and N ×1) forming a row and a column
of size N . In general, the inverse of a partitioned matrix B,
B =
(
P Q
R S
)
(A.14)
has the form [97]
B−1 =
(
P˜ Q˜
R˜ S˜
)
, (A.15)
115
where
P˜ = P−1 + (P−1 ·Q) · (S−R ·P−1 ·Q)−1 · (R ·P−1), (A.16)
Q˜ = −(P−1 ·Q) · (S−R ·P−1 ·Q)−1, (A.17)
R˜ = −(S−R ·P−1 ·Q)−1 · (R ·P−1), (A.18)
S˜ = (S−R ·P−1 ·Q)−1. (A.19)
These formulae can be easily checked by multiplying B with its inverse. The determi-
nant of a partitioned matrix can be computed through
detB = detP det(S−R ·P−1 ·Q) = detS det(P−Q · S−1 ·R). (A.20)
In our case, P = 1N , Q = (A1,N+1 . . . AN,N+1)T , R = (R1 . . . RN) and S = AN+1,N+1,
so the inverse takes the form
P˜ = 1+Q · (S−R ·Q)−1 ·R = 1+ 1λQ ·R
=


1 + λ−1A1,N+1R1 λ−1A1,N+1R2 · · · λ−1A1,N+1RN
λ−1A2,N+1R1 1 + λ−1A2,N+1R2 · · · λ−1A2,N+1RN
... ... . . . ...
λ−1AN,N+1R1 λ−1AN,N+1R2 · · · 1 + λ−1AN,N+1RN


,(A.21)
Q˜ = −Q · (S−R ·Q)−1 = −1λQ = −
1
λ(A1,N+1 . . . AN,N+1)
T , (A.22)
R˜ = −(S−R ·Q)−1 ·R = −1λR = −
1
λ(R1 . . . RN), (A.23)
S˜ = (S−R ·Q)−1 = 1λ, (A.24)
where λ ≡ det(S−R ·Q) = S−R ·Q = AN+1,N+1−
∑
i RiAi,N+1. From equation (A.20)
we see that λ = det(1 +∆A−1). This can also be obtained from the explicit expression
(A.13) by expanding around the last row and then the last column. Consequently, using
equation (A.11), λ is the ratio of determinants of A′ and A:
A′−1 = A−1(1 + ∆A−1)−1 ⇒ detA
′
detA =
detA−1
detA′−1 = det(1 + ∆A
−1) = λ (A.25)
Now it becomes clear why the inverse of the matrix A′ has to be computed – the terms
Ri in the expression for λ depend on all elements A−1k,i of the inverse of A. For an update
of A′ in turn we will need the inverse A′−1. Using equation (A.11) we finally arrive at
116
an explicit expression for A′−1:
A′−1 = A−1(1 + ∆A−1)−1 =


−λ−1L1
A−1i,j + λ−1LiRj
...
−λ−1LN
−λ−1R1 · · · −λ−1RN λ−1


,
(A.26)
where Li ≡
∑
k A−1i,kAk,N+1.
Let us now count the number of operations required for the calculation of detA′
and A′−1. For λ we need to compute the following quantities:
Ri =
∑
k AN+1,kA−1k,i : N2 mult. (N for each i), (N − 1)N add. (N − 1 for each i)∑
i RiAi,N+1 : N mult., N − 1 add.
So in total, to compute λ = AN+1,N+1 −
∑
i RiAi,N+1, N2 + N multiplications and N2
additions have to be performed. For the inverse matrix we need additionally:
Li =
∑
k A−1i,kAk,N+1 : N2 mult. (N for each i), (N − 1)N add. (N − 1 for each i)
λ−1 : 1 mult.
−λ−1Ri : N mult. (1 for each i)
−λ−1Li : N mult. (1 for each i)
A−1i,j + λ−1LiRj : N2 mult. (1 for each i and j), N2 add. (1 for each i and j)
In total this makes 2N2 + 2N + 1 multiplications and 2N2 −N additions and together
with the determinant 3N2 + 3N + 1 multiplications and 3N2 − N additions. So the
complexity of the algorithm is N2 and the prefactor is only 3 for both additions and
multiplications, which is better than the usual Sherman-Morrison algorithm that has a
prefactor of 4 in both cases.
A.3 Removing one row and one column
Assume we want to remove row n and column m from an N × N matrix A, which
inverse and determinant are already known. The resulting (N − 1) × (N − 1) matrix
A′ is extended to an N × N matrix in a similar manner as before, by replacing the
removed elements Anj by δmj and Aim by δin. Note that extending a matrix by a “unit”
row n and a “unit” column m implies extending the inverse matrix by a “unit” row m
and “unit” column n (this did not play a role for the previous update, where we had
117
m = n). Again as before (up to a minus sign) we define
∆ ≡ A− A′ =


...
0 Ai,m 0
...
· · · An,j · · · An,m − 1 · · · An,j · · ·
...
0 Ai,m 0
...


← n (A.27)
↑
m
and analogous to (A.11),
A′−1 = (A−∆)−1 = A−1(1−∆A−1)−1. (A.28)
Again we want to invert the term in brackets directly. Consider
∆A−1 =


...
0 Ai,m 0
...
· · · An,j · · · An,m − 1 · · · An,j · · ·
...
0 Ai,m 0
...


·


...
· · · A−1i,j · · ·...
· · · A−1m,j · · ·...
· · · A−1i,j · · ·...


=


...
· · · Ai,mA−1m,j · · ·...
· · · ∑k An,kA−1k,j − A−1m,j · · ·...
· · · Ai,mA−1m,j · · ·...


← n (A.29)
where the mth column of matrix ∆ and the mth row of matrix A−1 were highlighted, to
emphasise that the only contribution to product elements (∆A−1)ij in the rows i 6= n
comes from products Ai,mA−1m,j with i 6= n. The nth row of matrix ∆ (also highlighted)
118
multiplied with A−1 gives the nth row of the matrix ∆A−1 which has the form,
∑
k
An,kA−1k,j − A−1m,j = δnj − A−1m,j. (A.30)
Using this identity we can rewrite (A.29) as
∆A−1 =


...
· · · Ai,mA−1m,j · · ·...
· · · δnj − A−1m,j · · ·...
· · · Ai,mA−1m,j · · ·...


← n (A.31)
and consequently
1−∆A−1 =


...
· · · δij − Ai,mA−1m,j · · ·...
· · · A−1m,j · · ·...
· · · δij − Ai,mA−1m,j · · ·...


. ← n (A.32)
So the matrix 1 −∆A−1 consists of elements of the form δij − Ai,mA−1m,j, except row n
where the elements are of the form A−1m,j. The inverse of this expression, i.e. a matrix
M with ∑k(1 − ∆A−1)ikMkj = δij can be obtained by solving the following set of
equations:
∑
k
(δik − Ai,mA−1m,k)Mk,j = δij for i 6= n, (A.33)
∑
k
A−1m,kMk,j = δij for i = n. (A.34)
From (A.33) we obtain for i 6= n
Mi,j − Ai,m
∑
k
A−1m,kMk,j = δij (A.35)
119
and plugging in the identity ∑k A−1m,kMk,j = δnj from (A.34) we arrive at
Mi,j − Ai,mδnj = δij ⇔ Mi,j = δij + Ai,mδnj, (A.36)
for i 6= n. To determine the values of the remaining elements Mn,j we expand (A.34)
as follows,
δnj =
∑
k 6=n
A−1m,kMk,j + A−1m,nMn,j, (A.37)
and then replace Mk,j in the sum by the result from (A.36):
δnj =
∑
k 6=n
A−1m,k(δkj + Ak,mδnj) + A−1m,nMn,j, (A.38)
which we can do since the sum goes over all indices except n. It is now convenient to
consider two cases separately, j = n and j 6= n. For j = n,
δnn = 1 =
∑
k 6=n
A−1m,k(δkn + Ak,mδnn) + A−1m,nMn,n
=
∑
k 6=n
A−1m,k(0 + Ak,m) + A−1m,nMn,n
=
∑
k 6=n
A−1m,kAk,m + A−1m,nMn,n. (A.39)
Using
∑
k
A−1i,kAk,j = δij ⇒
∑
k 6=n
A−1m,kAk,m = δmm − A−1m,nAn,m = 1− A−1m,nAn,m (A.40)
equation (A.39) yields
1 = 1− A−1m,nAn,m + A−1m,nMn,n ⇔ Mn,n = An,m, (A.41)
provided A−1m,n 6= 0. However this is always the case for a non-singular matrix 1−∆A−1,
since if A−1m,n = 0, the elements of the nth column of 1−∆A−1 are all 0. For the second
case, j 6= n, (A.38) yields
0 =
∑
k 6=n
A−1m,k(δkj + Ak,m · 0) + A−1m,nMn,j = A−1m,j + A−1m,nMn,j. (A.42)
The sum could be performed as usual, since there exists a k 6= n such that δkj = 1 for
j 6= n. Altogether,
Mn,j = −
A−1m,j
A−1m,n
for j 6= n. (A.43)
120
Having calculated all elements of M ≡ (1 − ∆A−1)−1 we can write this matrix out
explicitly,
M = (1−∆A−1)−1 =


1 0 ...
. . . Ai,m 0
0 1 ...
· · · −A
−1
m,j
A−1m,n · · · An,m · · · −
A−1m,j
A−1m,n · · ·... 1 0
0 Ai,m . . .
... 0 1


← n(A.44)
↑
n
The determinant of M is the ratio of determinants of A and A′. It can be either obtained
by expandingM , first around the nth row and then each of the resulting matrices around
the rows that contain only one non-zero element, or by applying formula (A.20) for the
determinant of a partitioned matrix. Here the latter is demonstrated. We partition M
as follows: Let P be the n× n block in the upper left corner, S the (N − n)× (N − n)
unit matrix in the lower right corner and Q and R the n× (N −n) and the (N −n)×n
matrices in the upper right and lower left corner respectively. All elements of Q are 0
except the last row, with elements −A
−1
m,j
A−1m,n and similar for R where the non-zero elements
Ai,m are in the last column. So
detM = detS det(P−QS−1R) = det(P−Q ·R). (A.45)
The term Q ·R is an n× n matrix with only one non-zero element,
−
N∑
j=n+1
A−1m,j
A−1m,n
Aj,m,
in the lower right corner. So the matrix M˜ ≡ P−Q ·R has the form
M˜ =


...
1 Ai,m
...
· · · −A
−1
m,j
A−1m,n · · · An,m +
∑N
j=n+1
A−1m,j
A−1m,nAj,m


. (A.46)
121
To obtain the determinant of M˜ we apply (A.20) again,
det M˜ = det P˜ det(S˜− R˜P˜−1Q˜) = det(S˜− R˜ · Q˜), (A.47)
where the partitioning is analogous to the one used for the first update,
P˜ = 1n−1, (A.48)
S˜ = An,m +
N∑
j=n+1
A−1m,j
A−1m,n
Aj,m, (A.49)
R˜ = (· · · ,−A
−1
m,j
A−1m,n
, · · · ), (A.50)
Q˜ = (· · · , Ai,m, · · · )T . (A.51)
Altogether,
detM = det M˜ = An,m +
N∑
j=n+1
A−1m,j
A−1m,n
Aj,m +
n−1∑
j=1
A−1m,j
A−1m,n
Aj,m
= An,m +
∑
j 6=n
A−1m,jAj,m
A−1m,n
= An,m +
1− A−1m,nAn,m
A−1m,n
= 1A−1m,n
(A.52)
and thus
detA′
detA =
detA−1
detA′−1 =
1
detM = A
−1
m,n. (A.53)
What remains is to determine A′−1 = A−1(1−∆A−1)−1 by multiplying out the according
matrices. We obtain
A′−1i,j =
∑
k
A−1i,k (δkj −
A−1m,j
A−1m,n
δkn) = A−1i,j − A−1i,n
A−1m,j
A−1m,n
for j 6= n, (A.54)
A′−1i,j =
∑
k
A−1i,kAk,m = δim for j = n. (A.55)
Note that for i = m (A.54) simplifies to
A′−1m,j = A−1m,j − A−1m,n
A−1m,j
A−1m,n
= 0, (A.56)
122
so indeed A′−1 has “unit” row m and “unit” column n. These will be dropped at the
end of the procedure, yielding an updated matrix of size (N − 1)× (N − 1). Explicitly:
A′−1 =


...
A−1i,j − A−1i,n
A−1m,j
A−1m,n 0 A
−1
i,j − A−1i,n
A−1m,j
A−1m,n...
· · · 0 · · · 1 · · · 0 · · ·
...
A−1i,j − A−1i,n
A−1m,j
A−1m,n 0 A
−1
i,j − A−1i,n
A−1m,j
A−1m,n...


← m(A.57)
↑
n
Let us now count the number of operations for this update. There are no operations
required for the determinant. For each of the (N − 1)2 non-trivial elements of A′−1
we need to perform two multiplications and one addition. This can be performed even
faster, if we first calculate and store the N − 1 values of A
−1
m,j
A−1m,n for each j 6= n, which
reduces the number of required multiplications per element of A′−1 to one. The total
number of multiplications is thus N − 1 + (N − 1)2 = N2 −N and the total number of
additions is (N − 1)2 = N2 − 2N + 1, both of order N2 with prefactor 1. This update
is three times faster than the one that adds a row and a column.
A.4 Changing the elements of one row
The last update changes all elements in row n of A, for some n. The size of the matrix
remains unchanged. With the Sherman-Morrison scheme, this update requires order
2N2 additions and 2N2 multiplications. Here we apply a method, analogous to the ones
used for the previous updates and count the number of operations. First write
∆ ≡ A′ − A =


0 · · · 0
... ...
0 · · · 0
∆1 · · · ∆N
0 · · · 0
... ...
0 · · · 0


, ← n (A.58)
123
where ∆j ≡ A′n,j − An,j, and analogous to (A.11),
A′−1 = (A+∆)−1 = A−1(1 + ∆A−1)−1. (A.59)
It is easy to see that all entries of the matrix ∆A−1 are 0, except for row n where the
entries are
Rj ≡
∑
k
∆kA−1k,j =
∑
k
(A′n,k − An,k)A−1k,j =
∑
k
A′n,kA−1k,j − δnj. (A.60)
From this we get
1 + ∆A−1 =


1 0 0 · · · 0
. . . ... ...
0 1 0 · · · 0
· · · ∑k A′n,kA−1k,j · · ·
0 · · · 0 1 0
... ... . . .
0 · · · 0 0 1


← n (A.61)
Let us invert this matrix. For lucidity we introduce the abbreviation R′j ≡
∑
k A′n,kA−1k,j.
Since 1 + ∆A−1 varies from the unit matrix only in row n, then so does its inverse, or
in other words
(1 + ∆A−1)−1i,j = δij for i 6= n. (A.62)
To obtain the elements (1 + ∆A−1)−1n,j consider
δnj =
∑
k
R′k(1 + ∆A−1)−1k,j =
∑
k 6=n
R′kδk,j +R′n(1 + ∆A−1)−1n,j. (A.63)
Again there are two cases j = n and j 6= n. We deal with j = n first and obtain,
δnn = 1 =
∑
k 6=n
R′kδk,n +R′n(1 + ∆A−1)−1n,n = R′n(1 + ∆A−1)−1n,n. (A.64)
So, provided R′n 6= 0,
(1 + ∆A−1)−1n,n =
1
R′n
, (A.65)
but this is always the case for a non-singular matrix A′, since R′n = 0 would imply that
the nth column of 1 + ∆A−1 has only zero entries. For the case j 6= n:
δnj = 0 =
∑
k 6=n
R′kδk,j +R′n(1 + ∆A−1)−1n,j = R′j +R′n(1 + ∆A−1)−1n,j (A.66)
124
and hence
(1 + ∆A−1)−1n,j = −
R′j
R′n
. (A.67)
Explicitly:
(1 + ∆A−1)−1 =


1 0 0 · · · 0
. . . ... ...
0 1 0 · · · 0
· · · −R
′
j
R′n · · ·
1
R′n · · · −
R′j
R′n · · ·
0 · · · 0 1 0
... ... . . .
0 · · · 0 0 1


← n (A.68)
The determinant of this matrix is simply 1R′n as one can see by expanding around row
n (after removing row n, column n has only zero entries, so all minors will be zero
except the one where row n and column n are removed, in which case the minor is the
determinant of the unit matrix, i.e. 1). So the ratio of determinants of A′ and A is
detA′
detA =
detA−1
detA′−1 =
1
det(1 + ∆A−1)−1 = R
′
n. (A.69)
Now A′−1 can be obtained,
A′−1 = A−1(1 + ∆A−1)−1 =


...
A−1i,j − A−1i,n
R′j
R′n
A−1i,n
R′n A
−1
i,j − A−1i,n
R′j
R′n...

(A.70)
↑
n
Finally we need to count the number of operations. Each R′j ≡
∑
k A′n,kA−1k,j involves
N multiplications and N − 1 additions, which makes N2 multiplications and (N − 1)N
additions in total. Then we need to divide A−1i,n through R′n for each i. This requires
N multiplications. After this, for all i and j 6= n we need to multiply A
−1
i,n
R′n by R
′
j. The
total number of multiplications is thus N(N − 1). Finally we need N(N − 1) additions
to obtain the elements A−1i,j − A−1i,n
R′j
R′n . Altogether we have 2N
2 − 2N additions and
N2 + N + N(N − 1) = 2N2 multiplications. This is of order N2 with a prefactor
of 2 for both additions and multiplications, which is the same as we get using the
Sherman-Morrison method.
125
Bibliography
[1] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation. Oxford University
Press, 2003.
[2] C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases. Cam-
bridge University Press, 2002.
[3] S. Giorgini, L. P. Pitaevskii, and S. Stringari, “Theory of ultracold atomic Fermi
gases,” Rev. Mod. Phys. 80 no. 4, (Oct, 2008) 1215–1274, arXiv:0706.3360v2.
[4] I. Bloch, J. Dalibard, and W. Zwerger, “Many-body physics with ultracold gases,”
Rev. Mod. Phys. 80 no. 3, (Jul, 2008) 885–964, arXiv:0704.3011.
[5] O. Goulko and M. Wingate, “Monte Carlo study of a Fermi gas with infinite
scattering length,” PoS LAT2009 (2009) 062, arXiv:0910.3909.
[6] O. Goulko and M. Wingate, “Thermodynamics of balanced and slightly spin-
imbalanced Fermi gases at unitarity,” Phys. Rev. A 82 no. 5, (Nov, 2010) 053621,
arXiv:1008.3348.
[7] O. Goulko and M. Wingate, “The imbalanced Fermi gas at unitarity,” PoS
LAT2010 (2010) 187, arXiv:1011.0312.
[8] O. Goulko, F. Chevy, and C. Lobo, “Collision of two spin polarized fermionic
clouds,” to appear in Phys. Rev. A (Jun, 2011) , arXiv:1106.5773.
[9] J. J. Sakurai, Modern Quantum Mechanics. The Benjamin/Cummings Publishing
Company, Inc., 1985.
[10] L. D. Landau and E. M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory,
vol. 3 of Landau and Lifshitz: Course of Theoretical Physics. Pergamon Press,
1977.
[11] W. Ketterle and M. W. Zwierlein, “Making, probing and understanding ultracold
Fermi gases,” in Ultra-Cold Fermi Gases (Proceedings of the International School
of Physics “Enrico Fermi”, Course CLXIV), M. Inguscio, W. Ketterle, and C. Sa-
lomon, eds., pp. 95–288. IOS Press, 2007. arXiv:0801.2500.
126
[12] T. Ko¨hler, K. Go´ral, and P. S. Julienne, “Production of cold molecules via mag-
netically tunable Feshbach resonances,” Rev. Mod. Phys. 78 no. 4, (Dec, 2006)
1311–1361, cond-mat/0601420.
[13] S. Stringari, “Dynamics and superfluidity of an ultracold Fermi gas,” in Ultra-Cold
Fermi Gases (Proceedings of the International School of Physics “Enrico Fermi”,
Course CLXIV), M. Inguscio, W. Ketterle, and C. Salomon, eds., pp. 53–94. IOS
Press, 2007. cond-mat/0702526.
[14] D. M. Eagles, “Possible Pairing without Superconductivity at Low Carrier Con-
centrations in Bulk and Thin-Film Superconducting Semiconductors,” Phys. Rev.
186 (Oct, 1969) 456–463.
[15] P. Nikolic´ and S. Sachdev, “Renormalization-group fixed points, universal phase di-
agram, and 1/N expansion for quantum liquids with interactions near the unitarity
limit,” Phys. Rev. A 75 no. 3, (Mar, 2007) 033608, cond-mat/0609106.
[16] R. Haussmann, W. Rantner, S. Cerrito, and W. Zwerger, “Thermodynam-
ics of the BCS-BEC crossover,” Phys. Rev. A 75 no. 2, (Feb, 2007) 023610,
cond-mat/0608282.
[17] Y. Nishida and D. T. Son, “Fermi gas near unitarity around four and two spatial
dimensions,” Phys. Rev. A 75 no. 6, (Jun, 2007) 063617, cond-mat/0607835.
[18] Y. Nishida, “Unitary Fermi gas at finite temperature in the ² expansion,”
Phys. Rev. A 75 no. 6, (Jun, 2007) 063618, cond-mat/0608321.
[19] Y. Nishida and D. T. Son, “Unitary Fermi gas, epsilon expansion, and nonrelativis-
tic conformal field theories,” in The BCS-BEC crossover and the Unitary Fermi
Gas, W. Zwerger, ed., Lecture Notes in Physics. Springer-Verlag New York, LLC,
2011. arXiv:1004.3597.
[20] H. Hu, X.-J. Liu, and P. D. Drummond, “Temperature of a trapped uni-
tary Fermi gas at finite entropy,” Phys. Rev. A 73 no. 2, (Feb, 2006) 023617,
cond-mat/0509274.
[21] D. Lee and T. Scha¨fer, “Cold dilute neutron matter on the lattice. II. Results in
the unitary limit,” Phys. Rev. C 73 no. 1, (Jan, 2006) 015202, nucl-th/0509018.
[22] V. K. Akkineni, D. M. Ceperley, and N. Trivedi, “Pairing and superfluid properties
of dilute fermion gases at unitarity,” Phys. Rev. B 76 no. 16, (Oct, 2007) 165116,
cond-mat/0608154.
127
[23] A. Bulgac, J. E. Drut, and P. Magierski, “Quantum Monte Carlo simulations of
the BCS-BEC crossover at finite temperature,” Phys. Rev. A 78 no. 2, (Aug, 2008)
023625, arXiv:0803.3238.
[24] A. Bulgac, J. E. Drut, and P. Magierski, “Spin 1/2 Fermions in the Unitary Regime:
A Superfluid of a New Type,” Phys. Rev. Lett. 96 no. 9, (Mar, 2006) 090404,
cond-mat/0505374.
[25] A. Bulgac, J. E. Drut, and P. Magierski, “Thermodynamics of a Trapped Unitary
Fermi Gas,” Phys. Rev. Lett. 99 no. 12, (Sep, 2007) 120401, cond-mat/0701786.
[26] E. Burovski, N. Prokof’ev, B. Svistunov, and M. Troyer, “The Fermi-Hubbard
model at unitarity,” New J. Phys. 8 no. 8, (Aug, 2006) 153, cond-mat/0605350.
[27] L. Radzihovsky and D. E. Sheehy, “Imbalanced Feshbach-resonant Fermi gases,”
Rep. Prog. Phys. 73 no. 7, (Jul, 2010) 076501, arXiv:0911.1740.
[28] F. Chevy and C. Mora, “Ultra-cold polarized Fermi gases,” Rep. Prog. Phys. 73
no. 11, (Nov, 2010) 112401, arXiv:1003.0801.
[29] M. W. Zwierlein, A. Schirotzek, C. H. Schunck, and W. Ketterle, “Fermionic su-
perfluidity with imbalanced spin populations,” Science 311 no. 5760, (Jan, 2006)
492–496, cond-mat/0511197.
[30] S. Nascimbe`ne, N. Navon, K. J. Jiang, L. Tarruell, M. Teichmann, J. McKeever,
F. Chevy, and C. Salomon, “Collective Oscillations of an Imbalanced Fermi Gas:
Axial Compression Modes and Polaron Effective Mass,” Phys. Rev. Lett. 103
no. 17, (Oct, 2009) 170402, arXiv:0907.3032.
[31] G. B. Partridge, W. Li, R. I. Kamar, Y. an Liao, and R. G. Hulet, “Pairing and
Phase Separation in a Polarized Fermi Gas,” Science 311 no. 5760, (Jan, 2006)
503, cond-mat/0511752.
[32] M. M. Parish and D. A. Huse, “Evaporative depolarization and spin transport
in a unitary trapped Fermi gas,” Phys. Rev. A 80 no. 6, (Dec, 2009) 063605,
arXiv:0903.5324.
[33] X. Du, L. Luo, B. Clancy, and J. E. Thomas, “Observation of anomalous spin
segregation in a trapped Fermi gas,” Phys. Rev. Lett. 101 no. 15, (Oct, 2008)
150401, arXiv:0805.1036.
[34] X. Du, Y. Zhang, J. Petricka, and J. E. Thomas, “Controlling spin current
in a trapped Fermi gas,” Phys. Rev. Lett. 103 no. 1, (Jul, 2009) 010401,
arXiv:0901.3702.
128
[35] F. Pie´chon, J. N. Fuchs, and F. Laloe¨, “Cumulative identical spin rotation effects
in collisionless trapped atomic gases,” Phys. Rev. Lett. 102 no. 21, (May, 2009)
215301, arXiv:0901.4008.
[36] R. A. Duine, M. Polini, A. Raoux, H. T. C. Stoof, and G. Vignale, “Spin drag
in ultracold Fermi mixtures with repulsive interactions,” New J. Phys. 13 no. 4,
(Apr, 2011) 045010, arXiv:1012.3827.
[37] A. Sommer, M. Ku, G. Roati, and M. W. Zwierlein, “Universal Spin Trans-
port in a Strongly Interacting Fermi Gas,” Nature 472 (Apr, 2011) 201–204,
arXiv:1101.0780.
[38] G. M. Bruun, “Spin diffusion in Fermi gases,” New J. Phys. 13 no. 3, (Mar, 2011)
035005, arXiv:1012.1607.
[39] S. A. Wolf, D. D. Awschalom, R. A. Buhrman, J. M. Daughton, S. von Molna´r,
M. L. Roukes, A. Y. Chtchelkanova, and D. M. Treger, “Spintronics: A Spin-Based
Electronics Vision for the Future,” Science 294 no. 5546, (Nov, 2001) 1488–1495.
[40] A. Sommer, M. Ku, and M. W. Zwierlein, “Spin transport in polaronic and super-
fluid Fermi gases,” New J. Phys. 13 no. 5, (May, 2011) 055009, arXiv:1103.2337.
[41] J.-W. Chen and D. B. Kaplan, “Lattice theory for low energy fermions at
nonzero chemical potential,” Phys. Rev. Lett. 92 no. 25, (Jun, 2004) 257002,
hep-lat/0308016.
[42] D. Lee and R. Thomson, “Temperature-dependent errors in nuclear lattice simu-
lations,” Phys. Rev. C 75 no. 6, (Jun, 2007) 064003, nucl-th/0701048.
[43] E. Braaten, M. Kusunoki, and D. Zhang, “Scattering models for ultracold atoms,”
Annals Phys. 323 (Jul, 2008) 1770–1815, arXiv:0709.0499.
[44] M. E. Peskin and D. V. Schroeder, An Introduction to Quantum Field Theory.
Westview Press, 1995.
[45] E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics. Part 2, Theory of the
condensed state, vol. 9 of Landau and Lifshitz: Course of Theoretical Physics.
Pergamon Press, 1980.
[46] Y. Castin. Private communication.
[47] Y. Castin and F. Werner, “The Unitary Gas and its Symmetry Properties,” in The
BCS-BEC crossover and the Unitary Fermi Gas, W. Zwerger, ed., Lecture Notes
in Physics. Springer-Verlag New York, LLC, 2011. arXiv:1103.2851.
129
[48] Y. Castin, “Basic theory tools for degenerate Fermi gases,” in Ultra-Cold Fermi
Gases (Proceedings of the International School of Physics “Enrico Fermi”, Course
CLXIV), M. Inguscio, W. Ketterle, and C. Salomon, eds., pp. 289–349. IOS Press,
2007. cond-mat/0612613.
[49] E. M. Lifshitz and L. P. Pitaevskii, Physical Kinetics, vol. 10 of Landau and Lif-
shitz: Course of Theoretical Physics. Pergamon Press, 1981.
[50] A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, “Continuous-time quantum
Monte Carlo method for fermions,” Phys. Rev. B 72 no. 3, (Jul, 2005) 035122,
cond-mat/0411344.
[51] M. Campostrini, M. Hasenbusch, A. Pelissetto, and E. Vicari, “Theoretical es-
timates of the critical exponents of the superfluid transition in 4He by lattice
methods,” Phys. Rev. B 74 no. 14, (Oct, 2006) 144506.
[52] A. Pelissetto and E. Vicari, “Critical phenomena and renormalization-group the-
ory,” Phys. Rept. 368 no. 6, (Oct, 2002) 549 – 727, cond-mat/0012164.
[53] M. N. Barber, “Finite-size scaling,” in Phase Transitions and Critical Phenomena,
C. Domb and J. L. Lebowitz, eds., vol. 8, pp. 145–266. Academic Press, 1983.
[54] J. G. Brankov, Introduction to Finite-Size Scaling. Leuven University Press, 1996.
[55] D. Toussaint, “Simulating QCD at finite density,” Nucl. Phys. B (Proc. Suppl.)
17 (Sep, 1990) 248 – 251.
[56] S. Nascimbe`ne, N. Navon, S. Pilati, F. Chevy, S. Giorgini, A. Georges, and
C. Salomon, “Fermi-Liquid Behavior of the Normal Phase of a Strongly Inter-
acting Gas of Cold Atoms,” Phys. Rev. Lett. 106 no. 21, (May, 2011) 215303,
arXiv:1012.4664.
[57] Q. Chen, J. Stajic, S. Tan, and K. Levin, “BCS-BEC crossover: From high tem-
perature superconductors to ultracold superfluids,” Phys. Rep. 412 no. 1, (Jun,
2005) 1 – 88, cond-mat/0404274.
[58] T. DeGrand and C. DeTar., Lattice Methods for Quantum Chromodynamics. World
Scientific, 2006.
[59] S. Tan, “Energetics of a strongly correlated Fermi gas,” Ann. Phys. 323 no. 12,
(Dec, 2008) 2952 – 2970, cond-mat/0505200.
[60] S. Tan, “Large momentum part of a strongly correlated Fermi gas,” Ann. Phys.
323 no. 12, (Dec, 2008) 2971 – 2986, cond-mat/0508320.
130
[61] S. Tan, “Generalized virial theorem and pressure relation for a strongly correlated
Fermi gas,” Ann. Phys. 323 no. 12, (Dec, 2008) 2987 – 2990, arXiv:0803.0841.
[62] E. Braaten, “Universal Relations for Fermions with Large Scattering Length,” in
The BCS-BEC crossover and the Unitary Fermi Gas, W. Zwerger, ed., Lecture
Notes in Physics. Springer-Verlag New York, LLC, 2011. arXiv:1008.2922.
[63] F. Werner and Y. Castin, “Exact relations for quantum-mechanical few-body and
many-body problems with short-range interactions in two and three dimensions,”
(Jan, 2010) , arXiv:1001.0774.
[64] E. Burovski, E. Kozik, N. Prokof’ev, B. Svistunov, and M. Troyer, “Critical Tem-
perature Curve in BEC-BCS Crossover,” Phys. Rev. Lett. 101 no. 9, (Aug, 2008)
090402, arXiv:0805.3047.
[65] T. Abe and R. Seki, “From low-density neutron matter to the unitary limit,”
Phys. Rev. C 79 no. 5, (May, 2009) 054003, arXiv:0708.2524.
[66] T. Abe and R. Seki, “Lattice calculation of thermal properties of low-density neu-
tron matter with pionless NN effective field theory,” Phys. Rev. C 79 no. 5, (May,
2009) 054002, arXiv:0708.2523.
[67] H. Dong, L.-W. Siu, T. T. S. Kuo, and R. Machleidt, “Unitarity potentials and
neutron matter at the unitary limit,” Phys. Rev. C 81 no. 3, (Mar, 2010) 034003,
arXiv:0912.0273.
[68] J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, “Superfluid
Fermi Gases with Large Scattering Length,” Phys. Rev. Lett. 91 no. 5, (Jul, 2003)
050401, physics/0303094.
[69] G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, “Equation of State
of a Fermi Gas in the BEC-BCS Crossover: A Quantum Monte Carlo Study,”
Phys. Rev. Lett. 93 no. 20, (Nov, 2004) 200404, cond-mat/0406113.
[70] A. Gezerlis and J. Carlson, “Strongly paired fermions: Cold atoms and neutron
matter,” Phys. Rev. C 77 no. 3, (Mar, 2008) 032801, arXiv:0711.3006.
[71] J. T. Stewart, J. P. Gaebler, C. A. Regal, and D. S. Jin, “Potential Energy of a 40K
Fermi Gas in the BCS-BEC Crossover,” Phys. Rev. Lett. 97 no. 22, (Nov, 2006)
220406, cond-mat/0607776.
[72] L. Luo and J. E. Thomas, “Thermodynamic Measurements in a Strongly Interact-
ing Fermi Gas,” J. Low Temp. Phys. 154 (Jan, 2009) 1–29, arXiv:0811.1159.
131
[73] F. Palestini, A. Perali, P. Pieri, and G. C. Strinati, “Temperature and cou-
pling dependence of the universal contact intensity for an ultracold Fermi gas,”
Phys. Rev. A 82 no. 2, (Aug, 2010) 021605, arXiv:1005.1158.
[74] T. Enss, R. Haussmann, andW. Zwerger, “Viscosity and scale invariance in the uni-
tary Fermi gas,” Annals Phys. 326 no. 3, (Mar, 2011) 770 – 796, arXiv:1008.0007.
[75] H. Hu, X.-J. Liu, and P. D. Drummond, “Universal contact of strongly interacting
fermions at finite temperatures,” New J. Phys. 13 no. 3, (Mar, 2011) 035007,
arXiv:1011.3845.
[76] J. E. Drut, T. A. La¨hde, and T. Ten, “Momentum Distribution and Contact
of the Unitary Fermi Gas,” Phys. Rev. Lett. 106 no. 20, (May, 2011) 205302,
arXiv:1012.5474.
[77] S. Nascimbe`ne, N. Navon, K. J. Jiang, F. Chevy, and C. Salomon, “Exploring the
thermodynamics of a universal Fermi gas,” Nature 463 (Feb, 2010) 1057–1060,
arXiv:0911.0747.
[78] M. Horikoshi, S. Nakajima, M. Ueda, and T. Mukaiyama, “Measurement of Uni-
versal Thermodynamic Functions for a Unitary Fermi Gas,” Science 327 no. 5964,
(Jan, 2010) 442–445.
[79] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W. Zwierlein, “Revealing
the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary
Fermi Gas,” (Oct, 2011) , arXiv:1110.3309.
[80] Y. Shin, C. H. Schunck, A. Schirotzek, and W. Ketterle, “Phase diagram of a
two-component Fermi gas with resonant interactions,” Nature 451 (Feb, 2008)
689–693.
[81] N. Navon, S. Nascimbe`ne, F. Chevy, and C. Salomon, “The Equation of State of
a Low-Temperature Fermi Gas with Tunable Interactions,” Science 328 no. 5979,
(Apr, 2010) 729–732, arXiv:1004.1465.
[82] F. Werner. Private communication.
[83] X.-J. Liu, H. Hu, and P. D. Drummond, “Virial Expansion for a Strongly
Correlated Fermi Gas,” Phys. Rev. Lett. 102 no. 16, (Apr, 2009) 160401,
arXiv:0903.5366.
[84] K. Van Houcke, F. Werner, E. Kozik, N. Prokofev, B. Svistunov, M. Ku, A. Som-
mer, L. W. Cheuk, A. Schirotzek, and M. W. Zwierlein, “Feynman diagrams versus
Feynman quantum emulator,” (Oct, 2011) , arXiv:1110.3747.
132
[85] T.-L. Ho and E. J. Mueller, “High Temperature Expansion Applied to Fermions
near Feshbach Resonance,” Phys. Rev. Lett. 92 no. 16, (Apr, 2004) 160404,
cond-mat/0306187.
[86] S. Succi, The Lattice Boltzmann Equation. Oxford University Press, 2001.
[87] T. Lepers, D. Davesne, S. Chiacchiera, and M. Urban, “Numerical solution of the
Boltzmann equation for the collective modes of trapped Fermi gases,” Phys. Rev. A
82 no. 2, (Aug, 2010) 023609, arXiv:1004.5241.
[88] B. Jackson and E. Zaremba, “Modeling Bose-Einstein condensed gases at finite
temperatures with N -body simulations,” Phys. Rev. A 66 no. 3, (Sep, 2002)
033606, cond-mat/0205421.
[89] A. C. J. Wade, D. Baillie, and P. B. Blakie, “Direct simulation Monte Carlo method
for cold-atom dynamics: Classical Boltzmann equation in the quantum collision
regime,” Phys. Rev. A 84 no. 2, (Aug, 2011) 023612, arXiv:1105.2340.
[90] E. Taylor, S. Zhang, W. Schneider, and M. Randeria, “Colliding clouds of strongly
interacting spin-polarized fermions,” (Jun, 2011) , arXiv:1106.4245.
[91] S. Chiacchiera, T. Macr`ı, and A. Trombettoni, “Dipole oscillations in fermionic
mixtures,” Phys. Rev. A 81 no. 3, (Mar, 2010) 033624, arXiv:0911.4965.
[92] L. Vichi and S. Stringari, “Collective oscillations of an interacting trapped Fermi
gas,” Phys. Rev. A 60 no. 6, (Dec, 1999) 4734–4737, cond-mat/9905154.
[93] D. Gue´ry-Odelin, F. Zambelli, J. Dalibard, and S. Stringari, “Collective oscillations
of a classical gas confined in harmonic traps,” Phys. Rev. A 60 no. 6, (Dec, 1999)
4851–4856, cond-mat/9904409.
[94] M. Amoruso, I. Meccoli, A. Minguzzi, and M. P. Tosi, “Collective excitations of
a degenerate Fermi vapour in a magnetic trap,” Eur. Phys. J. D 7 no. 3, (Apr,
1999) 441–447, cond-mat/9907370.
[95] I. D’Amico and G. Vignale, “Theory of spin Coulomb drag in spin-polarized trans-
port,” Phys. Rev. B 62 no. 8, (Aug, 2000) 4853–4857.
[96] A. Trenkwalder, C. Kohstall, M. Zaccanti, D. Naik, A. I. Sidorov, F. Schreck,
and R. Grimm, “Hydrodynamic Expansion of a Strongly Interacting Fermi-Fermi
Mixture,” Phys. Rev. Lett. 106 no. 11, (Mar, 2011) 115304, arXiv:1011.5192.
[97] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical
Recipes in C++ : The Art of Scientific Computing. Cambridge University Press,
second ed., 2002.
133