DOCTORAL PROGRAM IN COMPUTATIONAL AND APPLIED PHYSICS UNIVERSITAT POLITÈCNICA
DE CATALUNYA DEPARTMENT OF PHYSICS Barcelona 2020

DOCTORAL THESIS
Author:  Supervisors:
Prof. Dr. Grigory Astrakharchik
Prof. Dr. Jordi Boronat Medico

Abstract In this Thesis, we report a detailed study of the ground-state
properties of a set of quantum few- and many-body systems by using Quantum
Monte Carlo methods. First, we introduced the Variational Monte Carlo and
Diffusion Monte Carlo methods, these are the methods used in this Thesis to
obtain the properties of the systems. The first systems we studied consist of
few-body clusters in a one-dimensional Bose-Bose and Bose-Fermi mixtures. Each
mixture is formed by two different species with attractive interspecies and
repulsive intraspecies contact interactions. For each mixture, we focused on
the study of the dimer, tetramer, and hexamer clusters. We calculated their
binding energies and unbinding thresholds. Combining these results with a
three-body theory, we extracted the three-dimer scattering length close to the
dimer-dimer zero crossing. For both mixtures, the three-dimer interaction
turns out to be repulsive. Our results constitute a concrete proposal for
obtaining a one-dimensional gas with a pure three-body repulsion. The next
system analyzed consists of few-body clusters in a two-dimensional Bose-Bose
mixture using two types of interactions. The first case corresponds to a
bilayer of dipoles aligned perpendicularly to the planes and, in the second,
we model the interactions by finite-range Gaussian potentials. We find that
all the considered clusters are bound states and that their energies are
universal functions of the scattering lengths, for sufficiently large
attraction-to-repulsion ratios. Studying the hexamer energy close to the
corresponding threshold, we discovered an effective three-dimer repulsion,
which can stabilize interesting many-body phases. Once the existence of the
bound states in the dipolar bilayer has been demonstrated, we investigated
whether halos can occur in this system. A halo state is a quantum bound state
whose size is much larger than the range of the attractive interaction between
the atoms that form it, showing universal ratios between energy and size. For
clusters composed from three up to six dipoles, we find two very distinct halo
structures. For large interlayer separation, the halo structure is roughly
symmetric. However, for the deepest bound clusters and as the clusters
approach the threshold, we discover an unusual shape of the halo states,
highly anisotropic. Importantly, our results prove the existence of stable
halo states composed of up to six particles. To the best of our knowledge,
this is the first time that halo states with such a large number of particles
have been predicted and observed in a numerical simulation. The next system we
studied is a two-dimensional many-body dipolar fluid confined to a bilayer
geometry. We calculated the ground-state phase diagram as a function of the
density and the separation between layers. Our simulations show that the
system undergoes a phase transition from a gas to a stable liquid as the
interlayer distance increases. The liquid phase is stable in a wide range of
densities and interlayer values. In the final part of this Thesis, we studied
a system of dipolar bosons confined to a multilayer geometry formed by equally
spaced two-dimensional layers. We calculated the ground-state phase diagram as
a function of the density, the separation between layers, and the number of
layers. The key result of our study in the dipolar multilayer is the existence
of three phases: atomic gas, solid, and gas of chains, in a wide range of the
system parameters. Remarkably, we find that the density of the solid phase
decreases several orders of magnitude as the number of layers in the system
increases. The results reported in this Thesis show that a dipolar system in a
bilayer and multilayer geometries offer stable and highly controllable setups
for observing interesting phases of quantum matter, such as halo states, and
ultra-dilute liquids and solids. Resumen En esta Tesis, presentamos un estudio
detallado de las propiedades del estado fundamental de un conjunto de sistemas
cuánticos de pocos y muchos cuerpos mediante el uso de los métodos de Monte
Carlo Cuántico. Primero, introducimos los métodos de Monte Carlo Variacional y
Monte Carlo Difusivo que usamos en esta Tesis para obtener las propiedades de
los sistemas. Los primeros sistemas que estudiamos son cúmulos de pocos
cuerpos en mezclas unidimensionales de Bose-Bose y Bose-Fermi. Cada una de las
mezclas está formada por dos especies con interacciones atractivas para
interespecies y repulsivas para intraespecies. Para cada una de las mezclas
nos enfocamos en el estudio de dímeros, tetrámeros y hexámeros. Calculamos las
energías de ligadura y los valores umbrales de ruptura de los cúmulos.
Combinando estos resultados con una teoría de tres cuerpos, extraemos la
longitud de dispersión de tres dímeros cerca del punto de ruptura dímero-
dímero. Para ambas mezclas la interacción de tres dímeros resulta ser
repulsiva. El siguiente sistema analizado son cúmulos de pocos cuerpos en una
mezcla bidimensional de Bose-Bose con dos tipos de interacciones. El primer
caso corresponde a una bicapa de dipolos con momentos dipolares orientados
perpendicularmente a los planos y, en el segundo, modelamos las interacciones
con potenciales gaussianos de rango finito. Encontramos que para relaciones de
atracción-repulsión suficientemente grandes todos los cúmulos considerados son
estados ligados y sus energías son funciones universales de las longitudes de
dispersión. Estudiando la energía del hexámero cerca del punto umbral
correspondiente, descubrimos una repulsión efectiva de tres dímeros, que puede
estabilizar fases interesantes de muchos cuerpos. Después de demostrar la
existencia de los estados ligados en la bicapa dipolar, investigamos si pueden
ocurrir estados de halo en este sistema. Un estado de halo es un estado ligado
cuántico cuyo tamaño es mucho mayor que el rango de la interacción atractiva
entre los átomos que lo forman. Para cúmulos compuestos de tres hasta seis
dipolos encontramos dos estructuras de halo muy distintas. Para separaciones
grandes entre las capas, la estructura de halo es aproximadamente simétrica.
Sin embargo, para los estados más ligados y a medida que los cúmulos se
acercan al punto umbral, descubrimos una estructura de halo inusual, altamente
anisotrópica. Nuestros resultados demuestran la existencia de estados de halo
estables compuestos de hasta seis dipolos. Hasta donde sabemos, esta es la
primera vez que estados de halo con un número tan grande de partículas se
predicen y observan en una simulación numérica. El siguiente sistema estudiado
es un fluido bidimensional dipolar de muchos cuerpos confinado a una geometría
de bicapa. Calculamos el diagrama de fases del estado fundamental como función
de la densidad y de la separación entre las capas. Nuestras simulaciones
muestran que en el sistema ocurre una transición de fase, de un gas a un
líquido a medida que se incrementa la distancia entre las capas. El líquido es
estable en una región amplia de densidades y de la distancia entre las capas.
En la parte final de esta Tesis se estudia un sistema de bosones dipolares
confinados a una geometría multicapa formada por capas bidimensionales
igualmente espaciadas. Calculamos el diagrama de fases del estado fundamental
como función de la densidad, la separación entre las capas y el número de
capas. El resultado clave de nuestro estudio sobre la multicapa es la
existencia de tres fases: gas atómico, sólido y gas de cadenas, en una región
amplia de los parámetros del sistema. Encontramos que la densidad del sólido
disminuye varios órdenes de magnitud a medida que el número de capas en el
sistema aumenta. Los resultados reportados en esta Tesis muestran que un
sistema de dipolos confinados a una bicapa o multicapa ofrecen configuraciones
estables y altamente controlables para observar fases interesantes de materia
cuántica.

###### Contents

  1. 1 Introduction
  2. 2 Quantum Monte Carlo Methods
    1. 2.1 The Quantum Many-Body Problem
    2. 2.2 Variational Monte Carlo Method
      1. 2.2.1 Variational Principle
      2. 2.2.2 The Method
      3. 2.2.3 Stochastic Processes
      4. 2.2.4 Markov Processes
      5. 2.2.5 Metropolis Algorithm
      6. 2.2.6 VMC Stochastic Realization
    3. 2.3 Diffusion Monte Carlo Method
      1. 2.3.1 Green’s Function
      2. 2.3.2 Importance Sampling
      3. 2.3.3 Importance-Sampling Green’s Function and Short-Time Approximation
      4. 2.3.4 DMC Stochastic Realization
      5. 2.3.5 Convergence Analysis
    4. 2.4 Trial Wave Functions
    5. 2.5 Quantum Monte Carlo Estimators
      1. 2.5.1 Pair Distribution Function
      2. 2.5.2 One-Body Density Matrix
      3. 2.5.3 Mixed Estimators and Extrapolation Technique
      4. 2.5.4 Pure Estimators
  3. 3 One-Dimensional Three-Boson Problem with Two- and Three-Body Interactions
    1. 3.1 Introduction
    2. 3.2 The System
    3. 3.3 Details of the Methods
    4. 3.4 Results
      1. 3.4.1 Binding Energies
      2. 3.4.2 Threshold Determination
      3. 3.4.3 Three-Dimer Repulsion
    5. 3.5 Summary
  4. 4 Few-Body Bound States of Two-Dimensional Bosons
    1. 4.1 Introduction
    2. 4.2 The Hamiltonian
    3. 4.3 Details of the Methods
    4. 4.4 Results
      1. 4.4.1 Binding Energies
      2. 4.4.2 Threshold Determination
      3. 4.4.3 Three-Dimer Repulsion
    5. 4.5 Summary
  5. 5 Quantum Halo States in Two-Dimensional Dipolar Clusters
    1. 5.1 Introduction
    2. 5.2 The Hamiltonian
    3. 5.3 Details of the Methods
    4. 5.4 Results
      1. 5.4.1 Structure of the Bound States
      2. 5.4.2 Quantum Halo Characteristics
    5. 5.5 Summary
  6. 6 Quantum Liquid of Two-Dimensional Dipolar Bosons
    1. 6.1 Introduction
    2. 6.2 The Hamiltonian
    3. 6.3 Details of the Methods
    4. 6.4 Results
      1. 6.4.1 Finite Size Effects
      2. 6.4.2 Equation of State
      3. 6.4.3 Phase Diagram
      4. 6.4.4 Depletion of the Condensate
      5. 6.4.5 Polarization
    5. 6.5 Summary
  7. 7 Phases of Dipolar Bosons Confined to a Multilayer Geometry
    1. 7.1 Introduction
    2. 7.2 The Hamiltonian
    3. 7.3 Details of the Methods
      1. 7.3.1 Gas Trial Wave Function
      2. 7.3.2 Solid Trial Wave Function
      3. 7.3.3 Gas of Chains Trial Wave Function
    4. 7.4 Details of the Gas of Chains Wave Function
      1. 7.4.1 Construction of Trial Wave Function
      2. 7.4.2 Many-Body Term: Gaussian Function
    5. 7.5 Results
      1. 7.5.1 Crystallization and Threshold Densities
      2. 7.5.2 Dipoles within a Four-Layer Geometry
      3. 7.5.3 Dipoles within an M-Layer Geometry
    6. 7.6 Summary
  8. 8 Conclusions
  9. A Trial Wave Function
  10. B Symmetric Trial Wave Function

###

##  Chapter 1 Introduction

In this Thesis, we report a detailed study of the ground-state properties of a
set of quantum few- and many-body systems in one and two dimensions with
different types of interactions. Nevertheless, the main focus of this work is
the study of the ground-state properties of an ultracold Bose system with
dipole-dipole interaction between the particles. We consider the cases where
the bosons are confined to a bilayer and multilayer geometries, that consist
of equally spaced two-dimensional layers. These layers can be experimentally
realized by imposing tight confinement in one direction. We specifically
address the study of new quantum phases, their properties, and transitions
between them. One expects these systems to have a rich collection of few- and
many-body phases because the dipole-dipole interaction is anisotropic and
quasi long-range. We will now present a short historical review of the
experiments and theoretical predictions that motivated the study of ultracold
dipolar bosonic gases.

The Bose-Einstein condensation (BEC) is a quantum phenomenon occurring when a
macroscopic number of bosons occupy the zero momentum state. This happens when
the system reaches a temperature below a critical value. Although BEC was
predicted by Albert Einstein in 1924 [1] based on a previous work by
Satyendranath Bose, it was not until 1995 that this phenomenon was
experimentally observed in rubidium [2] and sodium [3] gases independently.
The Nobel Prize in Physics 2001 was awarded to Wolfgang Ketterle, Eric A.
Cornell, and Carl E. Wieman for the achievement of BEC in dilute gases of
alkali atoms [4]. Typically the BEC state is reached for temperatures and
densities of the order of $T\sim 10^{-7}$K and $n\sim 10^{13}-10^{14}$cm-3,
respectively. Since the experimental observation of the BEC, there have been
intense theoretical and experimental efforts to understand ultracold bosonic
and fermionic gases. Interesting quantum phases have been predicted and
experimentally realized in these systems, for example, quantum droplets in a
mixture of Bose-Einstein condensates [5, 6, 7, 8] and in dipolar bosonic gases
[9, 10, 11, 12], quantum droplets in optical lattices [13, 14], impurities of
atoms immersed in a gas of fermions [15, 16, 17] or bosons [18, 19, 20], the
so-called polaron problem, among others. At very low densities, some ultracold
gases can be characterized by the s-wave scattering length, which means that
they can be described by an isotropic, short-range, contact interaction model.
However, there are gases with more complex interactions like dipolar
interactions.

Recent experiments have enabled the experimental study of ultracold gases with
dipole-dipole interaction (DDI). The DDI has two main properties that greatly
distinguish it from the contact interactions. Firstly, DDI is long-range in
three dimensions, it falls off with a power-law $1/r^{3}$ dependence, where
$r$ is the distance between particles. Secondly, DDI is anisotropic which
means that the interaction strength and its sign (repulsive or attractive),
depends on the angle between the polarization direction and the relative
distance of the particles. DDI can be found in magnetic atoms, ground-state
heteronuclear molecules, and Rydberg atoms, among others [21]. The first Bose-
Einstein condensate of magnetic atoms was realized in a gas of chromium atoms
in 2005 [22, 23]. The most recent experiments of dipolar gases are done with
Dysprosium [24, 25] and Erbium [26]. Many interesting phenomena have been
observed and predicted in dipolar gases, for example, dipolar Bose supersolid
stripes [27], dipolar quantum mixtures [28, 29], formation of a crystal phase
[30, 31], and a pair superfluid [32, 33, 34].

Describing a quantum many-body system is a demanding task, as it involves the
interactions of a large number of particles subject to spatial constraints.
Only for systems with very simple interactions, and under some assumptions,
can the Schrödinger equation be solved exactly. As we are studying systems
with dipolar interactions, complementary numerical methods become necessary,
like in our case, quantum Monte Carlo methods.

Quantum Monte Carlo (QMC) methods are a set of stochastic techniques that are
used to calculate the ground-state properties of quantum many-body systems at
zero or finite temperature [35, 36, 37]. One of the most used QMC techniques
for its simplicity is the Variational Monte Carlo (VMC) method. The VMC
technique uses the variational principle of quantum mechanics to provide an
upper bound to the ground-state energy of a quantum system. The accuracy of
this method depends entirely on the accuracy of the trial wave function used
to calculate the expectation value of the Hamiltonian. Another QMC technique
is the Diffusion Monte Carlo (DMC) method that solves the many-body
Schrödinger equation in imaginary time. This method consists of evolving in
imaginary time the wave function of a quantum system, and after enough time
has passed, it projects out the ground state. The DMC method allows one to
calculate the exact ground-state energy of the system, as well as other
properties, within controllable statistical errors. Both VMC and DMC methods
have been shown to give an accurate description of correlated quantum systems
[37]. Examples include ultradilute bosonic [38, 39] and fermionic mixtures
[29], Bose [40, 20] and Fermi [17] polarons, dipolar Bose supersolid stripes
[27], Bose gas subject to a multi-rod lattice [41], and ultracold quantum
gases with spin-orbit interactions [42].

In this Thesis, we have used QMC methods to study the ground-state properties
of a set of quantum few- and many-body systems. A large part of this Thesis is
focused on the study of dipolar Bose systems confined to a two-dimensional
bilayer and multilayer geometries. This Thesis is organized in the following
way:

In Chapter 2, we explain the basics of the Quantum Monte Carlo methods used in
this Thesis. First, we present the Variational Monte Carlo method, which is
used to calculate an approximation to the ground-state energy of a quantum
system. Then, we introduce the Metropolis algorithm, a method used to generate
random numbers from an arbitrary probability distribution function.
Afterwards, we discuss the Diffusion Monte Carlo method, which allows one to
calculate the exact ground-state energy of bosonic systems at zero
temperature. Later, we describe a number of trial wave functions used for QMC
calculations. Finally, we show how several ground-state properties are
evaluated in the Monte Carlo algorithm.

In Chapter 3, we use the DMC method to calculate the ground-state properties
of a one-dimensional Bose-Bose and Bose-Fermi mixtures with attractive
interspecies and repulsive intraspecies interactions. We focus on the study of
the tetramer and hexamer clusters. First, we describe the trial wave functions
for the system and the boundary conditions to be satisfied. Then, we evaluate
the tetramer and hexamer ground-state energies for Bose-Bose and Bose-Fermi
mixtures. Afterwards, we determine the threshold for unbinding for the
tetramer and hexamer, where the clusters break into two and three dimers,
respectively. Then, combining these results with a one-dimensional three-body
theory, we extract the three-dimer scattering length close to the dimer-dimer
zero crossing. Finally, we discuss a mixture of ultracold gases for obtaining
a one-dimensional gas with a pure three-body repulsion.

In Chapter 4, we study the ground-state properties of few-body bound states in
a two-dimensional mixture of A and B bosons with two types of interactions.
The first case corresponds to a bilayer of dipoles and, in the second, we
model the interactions by non-local (separable) finite-range Gaussian
potentials. First, we show the details of the numerical techniques used to
study the two models. In the dipolar case, we use the diffusion Monte Carlo
and in the Gaussian model, we use the stochastic variational method. Then,
using these methods we evaluate the ground-state binding energies of the
clusters. Also, we numerically determine the threshold for unbinding of the
bound states in the bilayer geometry. Afterwards, studying the hexamer energy
near to the tetramer threshold allows us to characterize an effective three-
dimer interaction, which may have important implications for the many-body
problem, particularly for observing liquid states of dipolar dimers in the
bilayer geometry. Finally, we give some examples of dipolar molecules as
promising candidates for observing the predicted few-body states within a
bilayer setup.

In Chapter 5, we analyze the ground-state properties of loosely bound dipolar
states confined to a two-dimensional bilayer geometry by using the VMC and DMC
methods. We study dipolar dimers, trimers, tetramers, pentamers, and hexamers.
First, we evaluate the pair distribution functions for the dimer, trimer, and
tetramer for different values of the interlayer separation. Then, we calculate
the spatial distributions functions for the trimer and tetramer for two
characteristic interlayer distances. Knowledge of these structural properties
permits us to understand how the size and shape of the clusters change with
the interlayer distance. Finally, the calculations of the binding energies and
sizes of the clusters allow us to investigate whether quantum halos, bound
states with a wave function that extends deeply into the classically forbidden
region, can occur in this system.

In Chapter 6, we study a many-body system of dipolar bosons within a bilayer
geometry by using exact many-body quantum Monte Carlo methods. We consider the
case in which the dipoles are aligned perpendicularly to the parallel layers.
First, we describe the trial wave functions for the system. Then, we calculate
the equation of state (energy per particle as a function of the density) for
different values of the interlayer distance. Knowledge of the equation of
state permits us to establish the quantum phases present in the bilayer of
dipoles. Afterwards, we obtain the gas-liquid phase diagram of the dipolar
fluid as a function of the density and the separation between layers. Finally,
we show numerical results for the one-body density matrix, condensate fraction
and polarization.

In Chapter 7, we use the diffusion Monte Carlo approach to study the ground-
state phase diagram of dipolar bosons in a geometry formed by equally spaced
two-dimensional layers. First, we discuss the trial wave functions to describe
the gas, solid, and gas of chains phases. In particular, for the trial
function of the chains, we have derived the expressions of the drift force and
the local energy, which are necessary to implement the DMC algorithm. Then, we
consider the case where there are four layers and the same number of dipoles
in each layer. In this case, we calculate the pair distribution functions for
the different phases present in the system. Also, we calculate the ground-
state phase diagram as a function of the total density and the interlayer
distance. Finally, we consider the case where the dipoles are confined to
three up to ten layers. Here, we calculate the zero-temperature phase diagram.

In Chapter 8, we present a summary of the principal results obtained in this
Thesis and the main conclusions achieved.

##### Pubications

The results of this doctoral research were published in:

  * •

G. Guijarro, A. Pricoupenko, G. E. Astrakharchik, J. Boronat, and D. S.
Petrov, One-dimensional three-boson problem with two- and three-body
interactions, Phys. Rev. A 97, 061605(R) (2018).

  * •

G. Guijarro, G. E. Astrakharchik, J. Boronat, B. Bazak, and D. S. Petrov, Few-
body bound states of two-dimensional bosons, Phys. Rev. A 101, 041602(R)
(2020).

Manuscripts in process

  * •

G. Guijarro, G. E. Astrakharchik, and J. Boronat, Quantum halo states in two-
dimensional dipolar clusters. Manuscript submitted for publication.

  * •

G. Guijarro, G. E. Astrakharchik, and J. Boronat, Quantum liquid of two-
dimensional dipolar bosons. Manuscript in preparation.

  * •

G. Guijarro, G. E. Astrakharchik, and J. Boronat, Phases of dipolar bosons
confined to a multilayer geometry. Manuscript in preparation.

##  Chapter 2 Quantum Monte Carlo Methods

The term Quantum Monte Carlo (QMC) refers to a set of stochastic techniques
whose objective is to solve as exactly as possible quantum many-body problems,
by determining the expectation values of quantum observables [35]. The QMC
methods have been demonstrated to give an accurate description of correlated
quantum systems at zero and low temperature [37]. Examples include ultracold
gases with bosonic [38, 39] and fermionic statistics [17, 29], quantum solids
[43, 44], and Helium [45, 46].

To study systems at zero temperature, one can use the Variational Monte Carlo
method (VMC) or the Diffusion Monte Carlo method (DMC). The VMC algorithm was
introduced by McMillan in 1965 to study liquid Helium [47]. In contrast, the
DMC technique was developed in several works over the years [48]. The VMC
method uses the variational principle of quantum mechanics to provide an upper
bound to the ground-state energy of a quantum system. On the other hand, the
DMC method allows one to calculate the _exact_ ground-state energy, of bosonic
systems by solving the many-body Schrödinger equation in imaginary time.

For fermionic systems, the DMC method provides an upper bound to the ground-
state energy and not the exact one [49]. This is because the wave function of
fermions is antisymmetric under the exchange of two particles. Therefore,
there are regions where it is positive and other regions where it is negative.
This leads to the so-called fermion sign problem.

To study quantum many-body systems with finite, but low temperature there
exists the Path Integral Monte Carlo (PIMC) method. This method is based on
the thermal density matrix and Feynman’s path-integral formulation of quantum
mechanics [50, 51].

In this chapter, we introduce the fundamental concepts of the Variational
Monte Carlo (VMC) and the Diffusion Monte Carlo (DMC) methods, which are the
Quantum Monte Carlo methods used in this Thesis. First, we discuss the
theoretical basis of the VMC method and its algorithm. Second, we present the
DMC method and its stochastic realization. Third, we discuss different types
of trial wave functions used in the method. Finally, we show how several
ground-state properties are evaluated in the Monte Carlo algorithm.

###  2.1 The Quantum Many-Body Problem

We will consider the generic quantum many-body problem involving $N$
interacting particles of mass $m$. We restrict ourselves to the case of
particles in an external potential $V_{\rm ext}(\mathbf{r}_{i})$ and pairwise
interactions $V_{\rm int}(\mathbf{r}_{i}-\mathbf{r}_{j})$. We can write the
Hamiltonian of such problems as

$\hat{H}=-\frac{\hbar^{2}}{2m}\sum_{i=1}^{N}\nabla^{2}_{\mathbf{r}_{i}}+\sum_{i=1}^{N}V_{\rm
ext}(\mathbf{r}_{i})+\sum_{i=1}^{N}\sum_{j=i+1}^{N}V_{\rm
int}(\mathbf{r}_{i}-\mathbf{r}_{j}),$ (2.1)

where $\mathbf{r}_{i}$ is the position of a single particle. It is difficult,
if not impossible, to exactly solve the Schrödinger equation for the many-body
Hamiltonian, which involves obtaining all its eigenstates. As the complete
analytical solution is unavailable, we use numerical methods to calculate the
wave function and the properties of the ground state. We would like to
calculate the ground-state expectation value of an observable $\hat{O}$

$\langle\hat{O}\rangle=\frac{\langle\Phi_{0}|\hat{O}|\Phi_{0}\rangle}{\langle\Phi_{0}|\Phi_{0}\rangle},$
(2.2)

$\Phi_{0}$ being the ground-state wave function. In particular, we are
interested in obtaining the ground-state energy of the system, which is
defined as

$E_{0}=\langle\hat{H}\rangle=\frac{\langle\Phi_{0}|\hat{H}|\Phi_{0}\rangle}{\langle\Phi_{0}|\Phi_{0}\rangle}.$
(2.3)

Using Monte Carlo methods we can calculate the exact value of the ground-state
energy of a Bose system at zero temperature, within some statistical errors.

###  2.2 Variational Monte Carlo Method

####  2.2.1 Variational Principle

The Variational Monte Carlo (VMC) method can be used to obtain an approximated
value of the ground-state energy of a quantum system by using the variational
principle of quantum mechanics. The variational principle states that the
expectation value of a Hamiltonian, $\hat{H}$, obtained with a trial wave
function $|\Psi_{\rm T}\rangle$, provides an upper bound to the ground-state
energy $E_{0}$ of the system:

$\frac{\langle\Psi_{\rm T}|\hat{H}|\Psi_{\rm T}\rangle}{\langle\Psi_{\rm
T}|\Psi_{\rm T}\rangle}\geq E_{0},$ (2.4)

if $|\Psi_{\rm T}\rangle$ is not orthogonal to the ground-state wave function.
The equality in Eq. (2.4) is fulfilled only when the trial function
$|\Psi_{\rm T}\rangle$ is the exact ground-state wave function. The proof of
Eq. (2.4) is as follows. If $|\phi_{n}\rangle$ is an eigenfunction with
eigenvalue $E_{n}$ of $\hat{H}$, the following properties are fulfilled

$\hat{H}|\phi_{n}\rangle=E_{n}|\phi_{n}\rangle,\qquad\langle\phi_{n}|\phi_{m}\rangle=\delta_{n,m},\quad{\rm
and}\quad\sum_{n}|\phi_{n}\rangle\langle\phi_{n}|=1.$ (2.5)

Using these relations the expectation value of $\hat{H}$ can be written as

$\displaystyle\frac{\langle\Psi|\hat{H}|\Psi\rangle}{\langle\Psi|\Psi\rangle}$
$\displaystyle=\frac{\sum_{n,m}\langle\Psi|\phi_{n}\rangle\langle\phi_{n}|\hat{H}|\phi_{m}\rangle\langle\phi_{m}|\Psi\rangle}{\sum_{n}\langle\Psi|\phi_{n}\rangle\langle\phi_{n}|\Psi\rangle}$
(2.6)
$\displaystyle=\frac{\sum_{n}E_{n}\langle\phi_{n}|\Psi\rangle|^{2}}{\sum_{n}|\langle\phi_{n}|\Psi\rangle|^{2}}.$

Since $E_{n}\geq E_{0}$, it follows that

$\displaystyle\frac{\langle\Psi|\hat{H}|\Psi\rangle}{\langle\Psi|\Psi\rangle}$
$\displaystyle=\frac{\sum_{n}E_{n}|\langle\phi_{n}|\Psi\rangle|^{2}}{\sum_{n}|\langle\phi_{n}|\Psi\rangle|^{2}}$
(2.7)
$\displaystyle\geq\frac{\sum_{n}E_{0}|\langle\phi_{n}|\Psi\rangle|^{2}}{\sum_{n}|\langle\phi_{n}|\Psi\rangle|^{2}}$
$\displaystyle=E_{0}\frac{\sum_{n}|\langle\phi_{n}|\Psi\rangle|^{2}}{\sum_{n}|\langle\phi_{n}|\Psi\rangle|^{2}}=E_{0},$

and this proves the upper bound reported in Eq. (2.4). In general, the trial
wave function $|\Psi_{\rm T}\rangle$ depends on a set of parameters that can
be optimized in order to find the lowest possible value of the energy. The
trial wave function with these optimal parameters is an approximation to the
ground-state wave function of $\hat{H}$ and the lowest energy is an upper
bound to the ground-state energy.

####  2.2.2 The Method

In the Variational Monte Carlo (VMC) method one defines a normalized
probability density function $\rho(\mathbf{R})$

$\rho(\mathbf{R})=\frac{|\Psi_{\rm T}(\mathbf{R})|^{2}}{\int
d\mathbf{R}|\Psi_{\rm T}(\mathbf{R})|^{2}},$ (2.8)

and a local energy $E_{\rm L}(\mathbf{R})$

$E_{\rm L}(\mathbf{R})=\frac{1}{\Psi_{\rm T}(\mathbf{R})}\hat{H}\Psi_{\rm
T}(\mathbf{R}),$ (2.9)

here $\mathbf{R}\leavevmode\nobreak\ =\leavevmode\nobreak\
(\vec{\mathbf{r}}_{1},\cdots,\vec{\mathbf{r}}_{N})$ is a 3$N$-dimensional
vector specifying the positions of $N$ particles. The expectation value of
$\hat{H}$ can be written in the integral form

$E_{\rm var}=\frac{\langle\Psi_{\rm T}|\hat{H}|\Psi_{\rm
T}\rangle}{\langle\Psi_{\rm T}|\Psi_{\rm T}\rangle}=\frac{\int
d\mathbf{R}\Psi_{\rm T}^{*}(\mathbf{R})\hat{H}\Psi_{\rm T}(\mathbf{R})}{\int
d\mathbf{R}|\Psi_{\rm T}(\mathbf{R})|^{2}}=\int
d\mathbf{R}\rho(\mathbf{R})E_{\rm L}(\mathbf{R}).$ (2.10)

The estimator of the variational energy $E_{\rm var}$ is then calculated as
the mean value of $E_{\rm L}(\mathbf{R})$:

$\overline{E_{\rm var}}=\frac{1}{M}\sum_{k=1}^{M}E_{\rm L}(\mathbf{R}_{k}),$
(2.11)

where $M$ is the number of points $\mathbf{R}_{k}$ sampled from the
probability density function $\rho(\mathbf{R})$. As we mentioned before,
$\Psi_{\rm T}$ depends on a set of parameters that are optimized to minimize
the energy. Therefore, we calculate the variational energy Eq. (2.11) for
several values of the parameters and obtain the minimum.

Other observables can also be calculated in the VMC method. The variational
expectation value of an observable $\hat{O}$ is given by

$\langle\hat{O}\rangle_{\rm var}=\frac{\int{\rm d}\mathbf{R}\Psi_{\rm
T}^{*}(\mathbf{R})\hat{O}\Psi_{\rm T}(\mathbf{R})}{\int{\rm
d}\mathbf{R}|\Psi_{\rm T}(\mathbf{R})|^{2}},$ (2.12)

which can be written as

$\langle\hat{O}\rangle_{\rm var}=\int{\rm d}\mathbf{R}\rho(\mathbf{R})O_{\rm
L}(\mathbf{R}),$ (2.13)

where $O_{\rm L}(\mathbf{R})$ is the local observable

$O_{\rm L}(\mathbf{R})=\frac{1}{\Psi_{\rm T}(\mathbf{R})}\hat{O}\Psi_{\rm
T}(\mathbf{R}).$ (2.14)

The variational estimator of any local observable can be computed by averaging
the corresponding local value

$\langle\hat{O}\rangle_{\rm var}=\frac{1}{M}\sum_{k=1}^{M}O_{\rm
L}(\mathbf{R}_{k}).$ (2.15)

In general, the probability density $\rho(\mathbf{R})=|\Psi_{\rm
T}(\mathbf{R})|^{2}/\int d\mathbf{R}|\Psi_{\rm T}(\mathbf{R})|^{2}$ Eq. (2.8)
is complicated and depends on many variables, thus it cannot be sampled by
using other methods such as the rejection method [48]. The solution to this
problem is found in the Metropolis algorithm which will be discussed below.
This method is used to generate random numbers from any probability
distribution by constructing a Markov process. Before presenting the
Metropolis algorithm, we are going to introduce the concepts of stochastic
processes and Markov processes.

####  2.2.3 Stochastic Processes

A stochastic process describes a time-dependent random variable
$\mathbf{R}(t)$. For times $t_{1},t_{2},\ldots,t_{n}$ there exist a
probability distribution

$P(\mathbf{R}_{1},t_{1};\mathbf{R}_{2},t_{2};\ldots;\mathbf{R}_{n},t_{n})$
(2.16)

where $\mathbf{R}_{1},\ldots,\mathbf{R}_{n}$ are random variables associated
to $\mathbf{R}(t)$. Usually the times are ordered, $t_{1}\leq
t_{2}\leq\ldots\leq t_{n}$. We can write the probability distribution in terms
of the conditional probabilities as

$\displaystyle
P(\mathbf{R}_{n},t_{n};\ldots;\mathbf{R}_{2},t_{2};\mathbf{R}_{1},t_{1})=$
$\displaystyle
P(\mathbf{R}_{n},t_{n}|\mathbf{R}_{n-1},t_{n-1};\ldots;\mathbf{R}_{1},t_{1})\ldots$
(2.17) $\displaystyle\times
P(\mathbf{R}_{2},t_{2}|\mathbf{R}_{1},t_{1})P(\mathbf{R}_{1},t_{1}).$

It is therefore clear that $\mathbf{R}_{j}$ is conditioned to
$\mathbf{R}_{j-1},\ldots,\mathbf{R}_{1}$. To calculate the probability
distribution of a particular realization of
$\mathbf{R}_{1},\ldots,\mathbf{R}_{n}$ we need to do it in order, this means,
first calculate $P(\mathbf{R}_{1},t_{1})$ then
$P(\mathbf{R}_{2},t_{2}|\mathbf{R}_{1},t_{1})$ and so on.

####  2.2.4 Markov Processes

A Markov process is a stochastic process for which the conditional probability
for the transition to a new state $\mathbf{R}_{j}$ depends only on the
previous state $\mathbf{R}_{j-1}$

$P(\mathbf{R}_{j},t_{j}|\mathbf{R}_{j-1},t_{j-1};\ldots;\mathbf{R}_{1},t_{1})=P(\mathbf{R}_{j},t_{j}|\mathbf{R}_{j-1},t_{j-1}).$
(2.18)

Therefore for a Markov process we can rewrite Eq. (2.17) as

$\displaystyle
P(\mathbf{R}_{n},t_{n};\ldots;\mathbf{R}_{2},t_{2};\mathbf{R}_{1},t_{1})=$
$\displaystyle P(\mathbf{R}_{n},t_{n}|\mathbf{R}_{n-1},t_{n-1})\ldots$ (2.19)
$\displaystyle\times
P(\mathbf{R}_{2},t_{2}|\mathbf{R}_{1},t_{1})P(\mathbf{R}_{1},t_{1}).$

From here onwards we will consider Markov processes independent of time which
are known as stationary Markov processes. The probability
$P(\mathbf{R}_{f},|\mathbf{R}_{i})$ is called the transition probability (or
matrix) of going from an initial state $\mathbf{R}_{i}$ to a final state
$\mathbf{R}_{f}$. The transition probability satisfy the following properties

$P(\mathbf{R}_{f}|\mathbf{R}_{i})\geq 0,$ (2.20) $\int
d\mathbf{R}_{f}P(\mathbf{R}_{f}|\mathbf{R}_{i})=1.$ (2.21)

The last property simply means that given an initial state $\mathbf{R}_{i}$, a
posterior state (the same or different) will be reached with certainty. Also,
there is not a fully absorbing state where the random walk stops.

We want to construct a Markov process that converges to the target probability
distribution $\rho(\mathbf{R})$ Eq. (2.8) by repeated applications of the
transition probability. In order for this to happen several conditions must be
met. The first one is that the distribution $\rho(\mathbf{R})$ must be an
eigenvector of $P(\mathbf{R}_{f}|\mathbf{R}_{i})$ with eigenvalue 1 [36]

$\int
d\mathbf{R}_{i}P(\mathbf{R}_{f}|\mathbf{R}_{i})\rho(\mathbf{R}_{i})=\rho(\mathbf{R}_{f})=\int
d\mathbf{R}_{i}P(\mathbf{R}_{i}|\mathbf{R}_{f})\rho(\mathbf{R}_{f})\quad\forall\mathbf{R}_{f},$
(2.22)

this condition is known as stationarity condition, which means that is we
start from the target distribution $\rho(\mathbf{R})$, after repeted
applications of the transition probability, we will continue to sample the
target distribution $\rho(\mathbf{R})$. In general it is required that
starting from any initial distribution $\rho_{\rm ini}(\mathbf{R})$, it should
converge to the target distribution $\rho(\mathbf{R})$ after applying the
transition probability a finite number of times,

$\displaystyle\lim_{n\to\infty}$ $\displaystyle\int
d\mathbf{R}_{1}d\mathbf{R}_{2}\ldots
d\mathbf{R}_{n}P(\mathbf{R}|\mathbf{R}_{n})P(\mathbf{R}_{n}|\mathbf{R}_{n-1})\ldots
P(\mathbf{R}_{2}|\mathbf{R}_{1})\rho_{\rm ini}(\mathbf{R}_{1})$ (2.23)
$\displaystyle=\rho(\mathbf{R}_{f}).$

To ensure the convergence to a unique stationary distribution
$\rho(\mathbf{R})$ the Markov process must be ergodic, which means that it
must be possible to move between any pair of states $\mathbf{R}_{j}$ and
$\mathbf{R}_{l}$ in a finite number of steps, then all the states can be
visited. Another condition that the Markov process must fulfill is the
detailed balanced contition

$P(\mathbf{R}_{f}|\mathbf{R}_{i})\rho(\mathbf{R}_{i})=P(\mathbf{R}_{i}|\mathbf{R}_{f})\rho(\mathbf{R}_{f}),$
(2.24)

for any states $\mathbf{R}_{i}$ and $\mathbf{R}_{f}$. This condition imposes
that the probability flux between the states $\mathbf{R}_{i}$ and
$\mathbf{R}_{f}$ to be the same in both directions.

####  2.2.5 Metropolis Algorithm

The Metropolis algorithm consists of a Markov process plus a decision
criterium on the random outcomes. We start with an initial state
$\mathbf{R}_{i}$. Then, we propose a temporary state $\mathbf{R}_{f}^{\prime}$
according to a probability distribution
$P_{\mathrm{prop}}(\mathbf{R}_{f}|\mathbf{R}_{i})$, which is known a priori.
After that, we test the temporary state. If the temporary state passes the
test then we accept it as the new initial state. If it does not pass the test
then the initial state remains unchanged. The test consists of accepting the
move (the temporary state) with probability
$P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})$ or rejecting the move with
probability $1-P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})$. Notice that,
the transition probability is given by

$P(\mathbf{R}_{f}|\mathbf{R}_{i})=\begin{cases}P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})P_{\mathrm{prop}}(\mathbf{R}_{f}|\mathbf{R}_{i})&\qquad\mathrm{if}\quad\mathbf{R}_{f}\neq\mathbf{R}_{i}\\\
1-\int
d\mathbf{R}^{\prime}_{f}P_{\mathrm{acc}}(\mathbf{R}^{\prime}_{f}|\mathbf{R}_{i})P_{\mathrm{prop}}(\mathbf{R}^{\prime}_{f}|\mathbf{R}_{i})&\qquad\mathrm{if}\quad\mathbf{R}_{f}=\mathbf{R}_{i}\end{cases}$
(2.25)

where $P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})$ is the probability of
accepting the move. We are free to choose the criterium for accepting a move,
this means we are free to choose
$P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})$. However,
$P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})$ has to fulfill the detailed
balanced condition Eq. (2.24)

$\frac{P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})}{P_{\mathrm{acc}}(\mathbf{R}_{i}|\mathbf{R}_{f})}=\frac{P_{\mathrm{prop}}(\mathbf{R}_{i}|\mathbf{R}_{f})\rho(\mathbf{R}_{f})}{P_{\mathrm{prop}}(\mathbf{R}_{f}|\mathbf{R}_{i})\rho(\mathbf{R}_{i})}.$
(2.26)

The Metropolis algorithm makes a particular choice of
$P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})$

$P_{\mathrm{acc}}(\mathbf{R}_{f}|\mathbf{R}_{i})=\mathrm{min}\left\\{1,\frac{P_{\mathrm{prop}}(\mathbf{R}_{i}|\mathbf{R}_{f})\rho(\mathbf{R}_{f})}{P_{\mathrm{prop}}(\mathbf{R}_{f}|\mathbf{R}_{i})\rho(\mathbf{R}_{i})}\right\\}.$
(2.27)

An advantage of this choice is that we do not need to calculate the
normalization factor for $\rho(\mathbf{R})$, because it will cancel out.

To implement the Metropolis algorithm we need to choose a proposal probability
$P_{\mathrm{prop}}(\mathbf{R}_{f}|\mathbf{R}_{i})$. A simple choice of
$P_{\mathrm{prop}}(\mathbf{R}_{f}|\mathbf{R}_{i})$ is a normal Gaussian
distribution.

The Metropolis algorithm reads as:

  1. 1.

Start from a random state $\mathbf{R}_{i}$.

  2. 2.

Propose a trial state $\mathbf{R}^{\prime}$ according to

$\mathbf{R}^{\prime}=\mathbf{R}_{i}+\mathbf{\chi},$

where $\mathbf{\chi}$ is an N-dimensional random vector sampled from a
Gaussian distribution.

  3. 3.

Calculate the quotient $\rho(\mathbf{R}^{\prime})/\rho(\mathbf{R}_{i})$.

  4. 4.

Generate a random number $\xi$ from the uniform distribution in $[0,1)$.

  5. 5.

If $\rho(\mathbf{R}^{\prime})/\rho(\mathbf{R}_{i})>\xi$ the move is accepted
and $\mathbf{R}_{i+1}=\mathbf{R}^{\prime}$. Otherwise stay in the same state
$\mathbf{R}_{i+1}=\mathbf{R}_{i}$.

After applying the Metropolis algorithm a large enough number of times, the
Markov process will sampled the target distribution $\rho(\mathbf{R})$.

Notice that in step 3 only the quotient
$\rho(\mathbf{R}^{\prime})/\rho(\mathbf{R}_{i})$ defines the acceptance
probabilty because
$P_{\mathrm{prop}}(\mathbf{R}_{i}|\mathbf{R}^{\prime})=P_{\mathrm{prop}}(\mathbf{R}^{\prime}|\mathbf{R}_{i})$,
since the Gaussian probability distribution is symmetric.

####  2.2.6 VMC Stochastic Realization

Here we present the VMC algorithm:

  1. 1.

We start with a random point $\mathbf{R}_{1}$ that represents the initial
distribution $\rho_{\rm ini}(\mathbf{R})$ given by

$\rho_{\rm ini}(\mathbf{R})=\delta(\mathbf{R}-\mathbf{R}_{1}).$

  2. 2.

Using the Metropolis algorithm we construct the Markov process given by
$\\{\mathbf{R}_{1},\mathbf{R}_{2},\ldots,\mathbf{R}_{B},\ldots,\mathbf{R}_{B+M}\\}$.

  3. 3.

We remove the first $B$ elements of the Markov process. The remaining elements
$\\{\mathbf{R}_{1},\mathbf{R}_{2},\ldots,\mathbf{R}_{M}\\}$ (with the
corresponding relabeling) are sampled from the target distribution
$\rho(\mathbf{R})$.

  4. 4.

Now we can calculate the variational estimator of the Hamiltonian

$\overline{E_{\rm var}}=\frac{1}{M}\sum_{k=1}^{M}E_{\rm L}(\mathbf{R}_{k}).$
(2.28)

###  2.3 Diffusion Monte Carlo Method

In the VMC method, the accuracy of the energy Eq. (2.28) depends entirely on
the accuracy of the trial wave function. The larger the overlap between the
trial wave function and the ground-state wave function the better the
estimation of the ground-state energy. To overcome the limitations of the VMC
method, we introduce the Diffusion Monte Carlo (DMC) method. This method
provides a practical way of evolving in imaginary time the wave function of a
quantum system and obtaining, ultimately, the ground-state energy [52].

The starting point of the DMC method is the time-dependent many-body
Schrödinger equation with an energy shift $E_{\rm T}$, which is equivalent to
replacing $\hat{H}\to\hat{H}-E_{\rm T}$

$\displaystyle i\hbar\frac{\partial\Psi(\mathbf{R},t)}{\partial t}$
$\displaystyle=\left[\hat{H}-E_{\rm T}\right]\Psi(\mathbf{R},t)$ (2.29)
$\displaystyle=\left[-\frac{\hbar^{2}}{2m}\nabla^{2}_{\mathbf{R}}+V(\mathbf{R})-E_{\rm
T}\right]\Psi(\mathbf{R},t),$

where $\mathbf{R}\leavevmode\nobreak\ =\leavevmode\nobreak\
(\vec{\mathbf{r}}_{1},\cdots,\vec{\mathbf{r}}_{N})$ is a 3$N$-dimensional
vector specifying the coordinates of all $N$ particles, $\Psi(\mathbf{R},t)$
is the many-body wave function of the system, which depends on the particle
coordinates and the time, and $\hat{H}$ is the many-body Hamiltonian Eq. (2.1)

$\nabla^{2}_{\mathbf{R}}=\sum_{i=1}^{N}\nabla^{2}_{\mathbf{r_{i}}},\quad
V(\mathbf{R})=\sum_{i=1}^{N}V_{\rm
ext}(\mathbf{r_{i}})+\sum_{i=1}^{N}\sum_{j=i+1}^{N}V_{\rm
int}(\mathbf{r_{i}}-\mathbf{r_{j}}).$ (2.30)

Let us now perform a transformation from real time to imaginary time by
introducing the new variable $\tau=it/\hbar$. After this, the Schrödinger
equation Eq. (2.29) becomes

$\displaystyle-\frac{\partial\Psi(\mathbf{R},\tau)}{\partial\tau}$
$\displaystyle=\left[\hat{H}-E_{\rm T}\right]\Psi(\mathbf{R},\tau)$ (2.31)
$\displaystyle=\left[-D\nabla^{2}_{\mathbf{R}}+V(\mathbf{R})-E_{\rm
T}\right]\Psi(\mathbf{R},\tau),$

where $D=\hbar^{2}/2m$. Eq. (2.31) can be identified as a modified diffusion
equation in the 3$N$\- dimensional space. If the $\left[V(\mathbf{R})-E_{\rm
T}\right]$ term were removed, Eq. (2.31) becomes the usual diffusion equation
with a difussion constant $D$. On the other hand, if the term with the
Laplacian were removed, Eq. (2.31) would be a rate equation, describing and
exponential growth or decrease of the function $\Psi(\mathbf{R},\tau)$.

The objective is to solve Eq. (2.31) to access the ground state of the system.
Using the spectral descomposition

$e^{-(\hat{H}-E_{\rm T})\tau}=\sum_{i}|\Phi_{i}\rangle e^{-(E_{i}-E_{\rm
T})\tau}\langle\Phi_{i}|,$ (2.32)

the formal solution of Eq. (2.31)

$|\Psi(\mathbf{R},\tau)\rangle=e^{-(\hat{H}-E_{\rm
T})\tau}|\Psi(\mathbf{R},0)\rangle,$ (2.33)

can be expressed as

$|\Psi(\mathbf{R},\tau)\rangle=\sum_{i=0}e^{-(E_{i}-E_{\rm
T})\tau}|\Phi_{i}\rangle\langle\Phi_{i}|\Psi(\mathbf{R},0)\rangle,$ (2.34)

where $\\{\Phi_{i}\\}$ and $\\{E_{i}\\}$, with
$\hat{H}|\Phi_{i}\rangle=E_{i}|\Phi_{i}\rangle$, denote a complete sets of
eigenfunctions and eigenvalues of $\hat{H}$, respectively. We consider that
the eigenvalues are ordered

$E_{0}<E_{1}\leq E_{2}\leq\ldots$ (2.35)

The amplitudes of each one of the terms in Eq. (2.34) can increase or decrease
in time depending on the sign of ($E_{n}-E_{\rm T}$). Notice that, for
sufficiently long times $\tau\to\infty$ the operator $e^{-(\hat{H}-E_{\rm
T})\tau}$ projects out the lowest eigenstate $|\Phi_{0}\rangle$ that has non-
zero overlap with $|\Psi(\mathbf{R},0)\rangle$

$\displaystyle\lim_{\tau\to\infty}|\Psi(\mathbf{R},\tau)\rangle=$
$\displaystyle\lim_{\tau\to\infty}\sum_{i=0}e^{-(E_{i}-E_{\rm
T})\tau}|\Phi_{i}\rangle\langle\Phi_{i}|\Psi(\mathbf{R},0)\rangle$ (2.36)
$\displaystyle=\lim_{\tau\to\infty}e^{-(E_{0}-E_{\rm
T})\tau}|\Phi_{0}\rangle\langle\Phi_{0}|\Psi(\mathbf{R},0)\rangle.$

The higher terms will decay exponentially faster since
$E_{n}>E_{0}\leavevmode\nobreak\ \leavevmode\nobreak\ \forall n\neq 0$. For
$E_{\rm T}=E_{0}$ the function $|\Psi(\mathbf{R},\tau)\rangle$ converges to
the ground-state wave function $|\Phi_{0}(\mathbf{R})\rangle$ regardless of
the choice of the initial wave function $|\Psi(\mathbf{R},0)\rangle$

$\lim_{\tau\to\infty}|\Psi(\mathbf{R},\tau)\rangle\propto|\Phi_{0}(\mathbf{R})\rangle.$
(2.37)

This fundamental property of the projector $e^{-(\hat{H}-E_{\rm T})\tau}$ is
the basis of the DMC technique [43]. The DMC method follows the evolution of
an initial many-body state $|\Psi(\mathbf{R},0)\rangle$ in imaginary time,
until long enough time passes and only the contribution of the ground state to
the many-body wave function dominates according to Eq. (2.36).

####  2.3.1 Green’s Function

To follow the evolution of the Schrödinger equation in imaginary time we will
use the Green’s function formalism.

The solution of the imaginary-time Schrödinger equation Eq. (2.31) in integral
form is given by

$\langle\mathbf{R}|\Psi(\tau)\rangle=\int
d\mathbf{R^{\prime}}\langle\mathbf{R}|e^{-(\hat{H}-E_{\rm
T})\tau}|\mathbf{R^{\prime}}\rangle\langle\mathbf{R^{\prime}}|\Psi(0)\rangle,$
(2.38)

and it can be written as

$\Psi(\mathbf{R},\tau)=\int
d\mathbf{R^{\prime}}G(\mathbf{R}|\mathbf{R^{\prime}};\tau)\Psi(\mathbf{R^{\prime}},0).$
(2.39)

Here, $\Psi(\mathbf{R^{\prime}},0)$ is the wave function at the initial time
$\tau=0$ and we have introduced the Green’s function
$G(\mathbf{R}|\mathbf{R^{\prime}};\tau)$, also known as the imaginary-time
propagator from $\mathbf{R^{\prime}}$ to $\mathbf{R}$

$G(\mathbf{R}|\mathbf{R^{\prime}};\tau)=\langle\mathbf{R}|e^{-(\hat{H}-E_{\rm
T})\tau}|\mathbf{R^{\prime}}\rangle.$ (2.40)

The Green’s function is subject to the boundary condition at the initial time
$\tau=0$

$G(\mathbf{R}|\mathbf{R^{\prime}};0)=\delta(\mathbf{R}-\mathbf{R^{\prime}}).$
(2.41)

In general, we do not know the exact Green’s function for all times $\tau$.
However, the Green’s function is known in the limit of a short propagation
time, $G(\mathbf{R}|\mathbf{R^{\prime}};\Delta\tau)$, where $\Delta\tau$ is a
small imaginary time-step

$\Psi(\mathbf{R},\tau+\Delta\tau)=\int
d\mathbf{R^{\prime}}G(\mathbf{R}|\mathbf{R^{\prime}};\Delta\tau)\Psi(\mathbf{R^{\prime}},\tau),$
(2.42)

and then Eq. (2.39) can be solved in a step by step process

$\displaystyle\Psi(\mathbf{R},\tau)=$ $\displaystyle\lim_{M\to\infty}\int
d\mathbf{R_{1}}d\mathbf{R_{2}}\cdots
d\mathbf{R_{M}}G(\mathbf{R}|\mathbf{R_{M}};\Delta\tau)G(\mathbf{R_{M}}|\mathbf{R_{M-1}};\Delta\tau)$
(2.43) $\displaystyle\cdots
G(\mathbf{R_{2}}|\mathbf{R_{1}};\Delta\tau)\Psi(\mathbf{R_{1}},0).$

According to Eq. (2.43) an approximation to the final state
$\Psi(\mathbf{R},\tau)$ is obtained by applying $M$ times the short-time
Green’s function to the intial state $\Psi(\mathbf{R_{1}},0)$.

Before giving an explicit expression for the short-time Green’s function we
are going to introduce the importance sampling technique. In this technique,
we introduce a guiding wave function that is independent of the imaginary
time.

####  2.3.2 Importance Sampling

Solving Eq. (2.29) is usually inefficient, mainly because of the presence of
the potential $V(\mathbf{R})$, which can diverge when to particles are very
close. This leads to large variance and low convergence when calculating the
expectation values of observables. To overcome these limitations one can use
the importance sampling technique.

In the importance sampling procedure one consider the imaginary-time evolution
of the mixed distribution $f(\mathbf{R},\tau)$, which is given by the product,

$f(\mathbf{R},\tau)=\Psi_{\rm T}(\mathbf{R})\Psi(\mathbf{R},\tau),$ (2.44)

of the wave function $\Psi(\mathbf{R},\tau)$, which satisfies the Schrödinger
equation Eq. (2.31), and a trial wave function $\Psi_{\rm T}(\mathbf{R})$,
which is imaginary-time independent. The trial wave function $\Psi_{\rm
T}(\mathbf{R})$ is designed from the available knowledge of the exact ground-
state wave function.

The imaginary-time evolution of $f(\mathbf{R},\tau)$ can be obtained by
multiplying Eq. (2.30) by $\Psi_{\rm T}(\mathbf{R})$. After rearranging terms,
one obtains

$\displaystyle-\frac{\partial f(\mathbf{R},\tau)}{\partial\tau}=$
$\displaystyle-D\nabla^{2}_{\mathbf{R}}f(\mathbf{R},\tau)+D\nabla_{\mathbf{R}}\cdot[\mathbf{F}(\mathbf{R})f(\mathbf{R},\tau)]$
(2.45) $\displaystyle+\left[E_{\rm L}(\mathbf{R})-E_{\rm
T}\right]f(\mathbf{R},\tau).$

Here, $\mathbf{F}(\mathbf{R})$ denotes the drift force, also called the drift
velocity

$\mathbf{F}(\mathbf{R})=2\frac{\nabla_{\mathbf{R}}\Psi_{\rm
T}(\mathbf{R})}{\Psi_{\rm T}(\mathbf{R})},$ (2.46)

and $E_{\rm L}(\mathbf{R})$ is the local energy Eq. (2.9)

$E_{\rm L}(\mathbf{R})=\frac{\hat{H}\Psi_{\rm T}(\mathbf{R})}{\Psi_{\rm
T}(\mathbf{R})}.$ (2.47)

Eq. (2.45) describes a modified difussion process for the mixed distribution
$f(\mathbf{R},\tau)$. Notice that, the rate term is now proportional to
$[E_{\rm L}(\mathbf{R})-E_{\rm T}]$, unlike the rate term in Eq. (2.31) which
depends on the potential $V(\mathbf{R})$. With a good choice of $\Psi_{\rm
T}(\mathbf{R})$, the local energy $E_{\rm L}(\mathbf{R})$ remain finite even
if $V(\mathbf{R})$ diverges [49]. Also, notice that there is an additional
term $\nabla_{\mathbf{R}}\cdot[\mathbf{F}(\mathbf{R})f(\mathbf{R},\tau)]$ in
Eq. (2.45). This new term imposes a drift on the difussion process guided by
$\Psi_{\rm T}(\mathbf{R})$.

The mixed distribution $f(\mathbf{R},\tau)$ becomes proportional to the
ground-state wave function in the limit of large enough time

$f(\mathbf{R},\tau)\propto{\lim_{\tau\to\infty}}\Psi_{\rm
T}(\mathbf{R})\Phi_{0}(\mathbf{R}).$ (2.48)

####  2.3.3 Importance-Sampling Green’s Function and Short-Time Approximation

The evolution described by Eq. (2.45) can be written as the sum of three
different operators acting on $f(\mathbf{R},\tau)$

$-\frac{\partial
f(\mathbf{R},\tau)}{\partial\tau}=(\hat{O}_{K}+\hat{O}_{D}+\hat{O}_{B})f(\mathbf{R},\tau)\equiv\hat{O}f(\mathbf{R},\tau),$
(2.49)

where

$\displaystyle\hat{O}_{K}$ $\displaystyle=-D\nabla^{2}_{\mathbf{R}},$ (2.50)
$\displaystyle\hat{O}_{D}$
$\displaystyle=D[\nabla_{\mathbf{R}}\cdot\mathbf{F}(\mathbf{R})+\mathbf{F}(\mathbf{R})\cdot\nabla_{\mathbf{R}}],$
$\displaystyle\hat{O}_{B}$ $\displaystyle=E_{\rm L}(\mathbf{R})-E_{\rm T}.$

Here, $\hat{O}_{K}$, $\hat{O}_{D}$ and $\hat{O}_{B}$ are the kinetic, drift
and branching operators, respectively. This division will make easier to solve
the Schrödinger equation for $f(\mathbf{R},\tau)$ Eq. (2.45).

Analogously to Eq. (2.42), the formal solution of the evolution equation for
the mixed distribution $f(\mathbf{R},\tau)$

$f(\mathbf{R},\tau+\Delta\tau)=\int
d\mathbf{R^{\prime}}\tilde{G}(\mathbf{R}|\mathbf{R^{\prime}};\Delta\tau)f(\mathbf{R^{\prime}},\tau),$
(2.51)

where $\tilde{G}(\mathbf{R}|\mathbf{R^{\prime}};\tau)$ is the importance
sampling Green’s function. $\tilde{G}(\mathbf{R}|\mathbf{R^{\prime}};\tau)$
satisfies the boundary condition

$\tilde{G}(\mathbf{R}|\mathbf{R^{\prime}};0)=\delta(\mathbf{R}-\mathbf{R^{\prime}}).$
(2.52)

The importance sampling Green’s function is given in terms of the operator
$\hat{O}$,

$\tilde{G}(\mathbf{R}|\mathbf{R^{\prime}};\Delta\tau)=\langle\mathbf{R}|e^{-\hat{O}\Delta\tau}|\mathbf{R^{\prime}}\rangle.$
(2.53)

Now, we focus on giving an explicit expression for the short-time Green’s
function. A short-time approximation of the Green’s function to first order in
$\Delta\tau$ is given by

$e^{-\hat{O}\Delta\tau}=e^{-(\hat{O}_{K}+\hat{O}_{D}+\hat{O}_{B})\Delta\tau}=e^{-\hat{O}_{K}\Delta\tau}e^{-\hat{O}_{D}\Delta\tau}e^{-\hat{O}_{B}\Delta\tau}+\mathcal{O}((\Delta\tau)^{2}).$
(2.54)

A second order descomposition is given by

$e^{-\hat{O}\Delta\tau}=e^{-\hat{O}_{B}\frac{\Delta\tau}{2}}e^{-\hat{O}_{D}\frac{\Delta\tau}{2}}e^{-\hat{O}_{K}\Delta\tau}e^{-\hat{O}_{D}\frac{\Delta\tau}{2}}e^{-\hat{O}_{B}\frac{\Delta\tau}{2}}+\mathcal{O}((\Delta\tau)^{3}).$
(2.55)

Observe that, as $\Delta\tau\rightarrow 0$ this will be a valid approximation.
Introducing Eq. (2.55) into Eq. (2.51) we obtain an integral equation of the
mixed distribution $f(\mathbf{R},\tau)$ in terms of the individual Green’s
functions $\tilde{G}_{i}$, each one asociated to a single operator
$\hat{O}_{i}$

$\displaystyle f(\mathbf{R},\tau+\Delta\tau)=$ $\displaystyle\int
d\mathbf{R_{1}}d\mathbf{R_{2}}d\mathbf{R_{3}}d\mathbf{R_{4}}d\mathbf{R^{\prime}}[\tilde{G}_{B}(\mathbf{R}|\mathbf{R_{1}};\frac{\Delta\tau}{2})$
(2.56)
$\displaystyle\times\tilde{G}_{D}(\mathbf{R_{1}}|\mathbf{R_{2}};\frac{\Delta\tau}{2})\tilde{G}_{K}(\mathbf{R_{2}}|\mathbf{R_{3}};\Delta\tau)$
$\displaystyle\times\tilde{G}_{D}(\mathbf{R_{3}}|\mathbf{R_{4}};\frac{\Delta\tau}{2})\tilde{G}_{B}(\mathbf{R_{4}}|\mathbf{R^{\prime}};\frac{\Delta\tau}{2})]f(\mathbf{R^{\prime}},\tau).$

The next step is to solve three diferential equations, each corresponding to a
Green’s function $\tilde{G}_{i}$. The first diferential equation is associated
with the kinetic operator $\tilde{G}_{K}$

$-\frac{\partial\tilde{G}_{K}(\mathbf{R}|\mathbf{R^{\prime}};\tau)}{\partial\tau}=-D\nabla^{2}_{\mathbf{R}}\tilde{G}_{K}(\mathbf{R}|\mathbf{R^{\prime}};\tau).$
(2.57)

This is a diffusion equation which diffusion constant $D$. The evolution given
by $\tilde{G}_{K}$ corresponds to an isotropic Gaussian movement

$\tilde{G}_{K}(\mathbf{R}|\mathbf{R^{\prime}};\tau)=\left(4\pi
D\tau\right)^{-3N/2}\mathrm{exp}\left[-\frac{(\mathbf{R}-\mathbf{R^{\prime}})^{2}}{4D\tau}\right].$
(2.58)

The second diferential equation corresponds to the drift operator
$\tilde{G}_{D}$

$-\frac{\partial\tilde{G}_{D}(\mathbf{R}|\mathbf{R^{\prime}};\tau)}{\partial\tau}=D\nabla_{\mathbf{R}}\cdot\left[\mathbf{F}(\mathbf{R})\tilde{G}_{D}(\mathbf{R}|\mathbf{R^{\prime}};\tau)\right].$
(2.59)

The Green’s function $\tilde{G}_{D}$ describes the movement due to the drift
force and its solution is

$\tilde{G}_{D}(\mathbf{R}|\mathbf{R^{\prime}};\tau)=\delta(\mathbf{R}-\mathcal{R}^{\prime}(\tau)),$
(2.60)

where $\mathcal{R}(\tau)$ is defined by the following equations

$\displaystyle\frac{d\mathcal{R}(\tau)}{d\tau}$
$\displaystyle=D\mathbf{F}(\mathcal{R}(\tau)),$ (2.61)
$\displaystyle\mathcal{R}(0)$ $\displaystyle=\mathbf{R}.$

The last diferential equation is associated with the branching operator
$\tilde{G}_{B}$

$-\frac{\partial\tilde{G}_{B}(\mathbf{R}|\mathbf{R^{\prime}};\tau)}{\partial\tau}=(E_{\rm
L}(\mathbf{R})-E_{\rm T})\tilde{G}_{B}(\mathbf{R}|\mathbf{R^{\prime}};\tau),$
(2.62)

and its solution is given by

$\tilde{G}_{B}(\mathbf{R}|\mathbf{R^{\prime}};\tau)=\mathrm{exp}\left[-(E_{\rm
L}(\mathbf{R})-E_{\rm
T})\tau\right]\delta\left(\mathbf{R}-\mathbf{R}^{\prime}\right).$ (2.63)

The Green’s function $\tilde{G}_{B}$ assigns a weight to $\mathbf{R}$
depending on its local energy.

Now that we have found the solutions to the equations of the Green’s functions
we can describe completely the stochastic realization of the DMC algorithm.

In the stochastic realization of the DMC algorithm, the mixed distribution and
its imaginary-time evolution are represented by a set of random walkers.
Walkers evolve through repeated applications of the propagators $G_{i}$, until
one obtains convergence to the ground state in the limit $\tau\to\infty$.

####  2.3.4 DMC Stochastic Realization

In this section, we use the concepts exposed previously to give a basic
version of the DMC algorithm with importance sampling.

In the DMC method, the probability distribution at the initial time
$f(\mathbf{R},0)$ and its evolution in imaginary-time $f(\mathbf{R},\tau)$ is
represented by a set of random walkers. A walker is defined by the positions
of all the particles in the system in the configuration space of 3$N$
dimensions
$\mathbf{R}=\\{\vec{\mathbf{r}}_{1},\vec{\mathbf{r}}_{2},\ldots,\vec{\mathbf{r}}_{N}\\}$.
The set of random walkers can be written as [41]

$\mathcal{R}_{k}=\\{\mathbf{R}_{k,\alpha}|\alpha=1,2,\ldots,N_{w,k}\\}.$
(2.64)

Here, $k$ is the time step index, $\tau_{k}=k\Delta\tau$ is the current time,
and $N_{w,k}$ is the number of walkers which may change between steps.

The initial configuration for the DMC algorithm is drawn from some arbitrary
probability distribution. In most cases the initial configuration will be the
output from the VMC algorithm.

At the time-step $k=0$, we start with an initial configuration $N_{w,0}$ of
random walkers $\mathcal{R}_{0}$ drawn from

$f_{\rm ini}(\mathbf{R})=f(\mathbf{R},0)=|\Psi_{\rm T}(\mathbf{R})|^{2},$
(2.65)

after passing a large enough Metropolis steps. An initial estimate of $E_{\rm
T}$ is obtained from the mean of the local energies of the walkers

$E_{T,0}=\frac{1}{N_{w,0}}\sum_{\alpha=1}^{N_{w,0}}E_{\rm
L}(\mathbf{R}_{0,\alpha}).$ (2.66)

The following algorithm is iterated $M$ times.

  1. 1.

Starting from the random walker $\mathbf{R}_{k-1,\alpha}$ we obtain a
temporary configuration by appliyng the Green’s function $\hat{G}_{K}$ Eq.
(2.57). This is done for all random walkers in $\mathcal{R}_{k-1}$. This means
that if we start with the configuration $\mathbf{R}_{k-1,\alpha}$ we obtain a
temporary configuration $\mathbf{R}^{\prime}_{k-1,\alpha}$ as

$\mathbf{R}^{\prime}_{k-1,\alpha}=\mathbf{R}_{k-1,\alpha}+\mathbf{\chi}.$
(2.67)

Here, $\mathbf{\chi}$ is an N-dimensional random vector sampled from a
multivariate Gaussian distribution with zero mean and variance
$\sigma^{2}=2D\Delta\tau$.

  2. 2.

Now we will apply a second Green’s function $\hat{G}_{D}$ Eq. (2.59), which
corresponds to the action of the drift force. From the temporary configuration
$\mathbf{R^{\prime}}_{k-1,\alpha}$ we obtain a new configuration
$\mathbf{R}_{k,\alpha}$ by doing the following steps

     * •

$\mathbf{R}^{(1)}_{k-1,\alpha}=\mathbf{R}^{\prime}_{k-1,\alpha}+\mathbf{F}(\mathbf{R}^{\prime}_{k-1,\alpha})\frac{\Delta\tau}{2}$

     * •

$\mathbf{R}^{(2)}_{k-1,\alpha}=\mathbf{R}^{\prime}_{k-1,\alpha}+\left[\mathbf{F}(\mathbf{R}^{\prime}_{k-1,\alpha})+\mathbf{F}(\mathbf{R}^{(1)}_{k-1,\alpha})\right]\frac{\Delta\tau}{4}$

     * •

$\mathbf{R}^{(3)}_{k-1,\alpha}=\mathbf{R}^{\prime}_{k-1,\alpha}+\mathbf{F}(\mathbf{R}^{(2)}_{k-1,\alpha})\Delta\tau$

     * •

$\mathbf{R}_{k,\alpha}=\mathbf{R}^{(3)}_{k-1,\alpha}$

The new configuration form a new set $\mathcal{R}_{k}$.

We used a second order integration method called Runge-Kutta [53] to
integrated the differential equation Eq. (2.61) in order to do the
displacement from $\mathbf{R^{\prime}}_{k-1,\alpha}$ to
$\mathbf{R}_{k,\alpha}$.

  3. 3.

Calculate the branching probability for each walker in $\mathcal{R}_{k}$:

$w_{\alpha}=e^{-\left(\frac{E_{\rm L}(\mathbf{R}_{k,\alpha})-E_{\rm
L}(\mathbf{R}_{k-1,\alpha})}{2}-E_{\rm T}\right)\Delta\tau}.$ (2.68)

  4. 4.

Calculate the branching factor for each walker in $\mathcal{R}_{k}$:

$n_{\alpha}=int(w_{\alpha}+\eta).$ (2.69)

Here, $int$ denotes the integer part of a real number and $\eta$ is a random
number drawn from the uniform distribution on the interval $[0,1)$. If
$n_{\alpha}=0$, removed $\mathbf{R}_{k,\alpha}$ from $\mathcal{R}_{k}$. If
$n_{\alpha}\geq 1$, replace $\mathbf{R}_{k,\alpha}$ with $n_{\alpha}$ copies
of iftsel in $\mathcal{R}_{k}$.

  5. 5.

Update the estimators of energy and other observables of interest.

  6. 6.

Repeat steps 1 to 5 until $M$ time steps are reached.

For sufficiently long times the ground-state energy is given by

$E_{0}=\lim_{\tau\to\infty}\frac{\int{d}\mathbf{R}f(\mathbf{R},\tau)E_{\rm
L}(\mathbf{R})}{\int d\mathbf{R}f(\mathbf{R},\tau)}.$ (2.70)

The result of the stochastic process describe above is a set of walkers
representing the distribution $f(\mathbf{R},\tau)$. Therefore, the estimator
for $\overline{E_{0}}$ after $M$ times steps is [36]

$\overline{E_{0}}=\frac{\sum_{k=1}^{M}\sum_{\alpha=1}^{N_{w,k}}E_{\rm
L}(\mathbf{R}_{k,\alpha})}{\sum_{k=1}^{M}\sum_{\alpha=1}^{N_{w,k}}w_{k,\alpha}}.$
(2.71)

The value of $E_{\rm T}$ is adjusted during the iterations to keep the size of
the walker population within a desired value. A simple formula for adjusting
$E_{\rm T}$ for the iteration $k+1$ is [36]

$E_{T,k+1}=\overline{E_{0,k}}-C\mathrm{ln}\left(\frac{N_{w,k}}{N_{w,ave}}\right),$
(2.72)

where $C$ is a constant, and $N_{k,ave}$ is a desired average number of
walkers.

####  2.3.5 Convergence Analysis

The DMC algorithm gives exact results for the ground-state energy when
simultaneously the time step $\Delta\tau\to 0$ and the number of walkers
$N_{w}\to\infty$. The use of a finite time step $\Delta\tau$ to approximate
the Green’s function introduces a systematic error bias in the calculation. To
overcome this problem one can consider a short-time Green’s function accurate
to order $(\Delta\tau)^{2}$ according to Eq. (2.54) or a more precise
algorithm accurate to order $(\Delta\tau)^{3}$ as stated by Eq. (2.55). In the
first case, the energy has a linear dependence when the time step is
sufficiently small. Then, one can use several values of the time step to
extrapolate the value of the energy to the $\Delta\tau\to 0$ limit. In the
second case, the energy depends quadratically on the time step. Here, the
extrapolation procedure is not completely necessary because for small
$\Delta\tau$ the energy converges fast to the exact value and the time step
can be chosen such that the systematic error is smaller than the statistical
error.

In Fig 2.1 we show an example of the ground-state energy $E_{0}$ dependence on
the time step $\Delta\tau$ for a dipolar gas. We compare a linear DMC method
with a second-order DMC method. For the linear DMC algorithm, in order to
obtain the exact energy, the extrapolation to zero time step is required. In
Fig 2.1 we observe that the slope of the line that joins the green dots is
pronounced with respect to the scale that we are using. In contrast, for the
second-order DMC technique, we notice that the changes in energy are
statistically indistinguishable in the range $\Delta\tau=0.01-0.1$, the slope
of the line that joins the blue dots is less pronounced. This is a very useful
feature of the second-order DMC method which allows one to obtain exact values
of the energy with less computational effort.

Figure 2.1: Ground-state energy $E_{0}$ of a dipolar gas vs the time step
$\Delta\tau$ computed by using the linear and second-order DMC algorithms. The
dashed lines correspond to polynomial fits:
$E_{0}(\Delta\tau)=a_{0}+a_{1}\Delta\tau$ for linear DMC meyhod and
$E_{0}(\Delta\tau)=b_{0}+b_{1}\Delta\tau+b_{2}(\Delta\tau)^{2}$ for second-
order DMC method. The unit of energy used is $\hbar^{2}/(mr_{0}^{2})$, with
$r_{0}=md^{2}/\hbar^{2}$ the dipolar length and $d$ is the dipole moment of an
atom of mass $m$. The units of $\Delta\tau$ are inverse of energy.

Besides the time step bias, the DMC algorithm presents a dependence on the
number of walkers $N_{w}$, which requires additional analysis.

Fig 2.2 shows an example of the ground-state energy $E_{0}$ dependence on the
number of walkers $N_{w}$ for a dipolar gas. We notice that the energy changes
very little for $N_{w}\gtrsim 1000$. Therefore, using $N_{w}\approx 1000$ to
estimate the exact ground-state energy is typically sufficient. This depends
on the interaction potential and mainly on the quality of the trial wave
function. The improvement of $\Psi_{\rm T}$ makes $N_{w}$ decrease.

Figure 2.2: Ground-state energy $E_{0}$ of a dipolar gas vs the the inverse of
the number of walkers $1/N_{w}$ computed by using a second-order DMC algorithm
with time step $\Delta\tau=0.01$. The dashed line correspond to the polynomial
fit: $E_{0}(N_{w})=a_{0}+a_{1}/N_{w}+a_{2}/(N_{w})^{2}+a_{3}/(N_{w})^{3}.$ The
unit of energy used is $\hbar^{2}/(mr_{0}^{2})$. The units of $\Delta\tau$ are
inverse of energy.

###  2.4 Trial Wave Functions

An important part of the VMC and DMC methods is the choice of the trial wave
function. In the VMC technique, the expectation values of all observables are
evaluated with the trial wave function, therefore it determines completely the
accuracy of the results. In the DMC algorithm, the trial function affects the
efficiency of the estimations by increasing or decreasing the variance. The
DMC technique is based on energy projection, therefore, generally, the energy
converges faster as compared to other quantities. However, the non-diagonal
properties can be more sensitive to the quality of the trial wave function.

The trial wave function $\Psi_{\rm T}(\mathbf{R})$ should be a good
approximation of the ground state of the system. Also, for better
computational efficiency $\Psi_{\rm T}(\mathbf{R})$ and its gradient and
Laplacian should have simple expressions since they are repeatedly evaluated
in the calculation.

The trial wave functions usually used in Quantum Monte Carlo methods are of
the form

$\Psi_{\rm
T}(\mathbf{R})=\mathcal{F}_{1}(\mathbf{R})\mathcal{F}_{2}(\mathbf{R})\mathcal{S}(\mathbf{R}).$
(2.73)

The factor $\mathcal{F}_{1}(\mathbf{R})$ is constructed as a product of one-
body terms

$\mathcal{F}_{1}(\mathbf{R})=\prod_{i=1}^{N}f_{1}(\mathbf{r_{i}}).$ (2.74)

The one-body term $f_{1}(\mathbf{r_{i}})$ depends only on the position of a
single particle $\mathbf{r_{i}}$. The choice of $f_{1}(\mathbf{r_{i}})$ is
based on the characteristics of the system under study. In general, the one
body functions are taken as the solution of the one-body problem with an
external potential $V_{ext}(\mathbf{r_{i}})$.

The interparticle correlations are commonly described by a pair-product form
$\mathcal{F}_{2}(\mathbf{R})$ known as the Bijl-Jastrow term and it is
constructed as a product of two-body terms

$\mathcal{F}_{2}(\mathbf{R})=\prod_{j<k}^{N}f_{2}(|\mathbf{r_{j}}-\mathbf{r_{k}}|).$
(2.75)

The two-body term depends on the distance between a pair of particles
$f_{2}(|\mathbf{r_{j}}-\mathbf{r_{k}}|)$. Typically, for short distances the
two-body function is constructed as the solution of the two-body problem with
an interaction potential $V_{int}(|\mathbf{r_{j}}-\mathbf{r_{k}}|)$.

The factor $\mathcal{S}(\mathbf{R})$ defines the symmetry or antisymmetry of
$\Psi_{\rm T}(\mathbf{R})$ under the exchange of two particles.

In Chapters 4, 5, and 6 we study a mixture of bosons of types A and B with
attractive interspecies AB interactions and equally repulsive intraspecies AA
and BB interactions. In this case, we use $\Psi_{\rm S}(\mathbf{R})$ as a
trial wave function

$\displaystyle\Psi_{\rm S}(\mathbf{R})=$
$\displaystyle\prod_{i<j}^{N_{\text{A}}}f_{{\text{A}}{\text{A}}}(|\mathbf{r}_{i}-\mathbf{r}_{j}|)\prod_{\alpha<\beta}^{N_{{\text{B}}}}f_{{\text{B}}{\text{B}}}(|\mathbf{r}_{\alpha}-\mathbf{r}_{\beta}|)$
(2.76)
$\displaystyle\times\left[\prod_{i=1}^{N_{{\text{A}}}}\sum_{\alpha=1}^{N_{{\text{B}}}}f_{{\text{A}}{\text{B}}}(|\mathbf{r}_{i}-\mathbf{r}_{\alpha}|)+\prod_{\alpha=1}^{N_{{\text{B}}}}\sum_{i=1}^{N_{{\text{A}}}}f_{{\text{A}}{\text{B}}}(|\mathbf{r}_{i}-\mathbf{r}_{\alpha}|)\right],$

where $N_{\text{A}}$ and $N_{\text{B}}$ are the number of bosons of the
species A and B, respectively. We denote with Latin letters the bosons of the
species A and with Greek letters the bosons of the species B. In Eq. (2.76) we
have removed the one-body terms since there is no a external potential. The
terms in the first row of Eq. (2.76) are of the Bijl-Jastrow form and the term
in the second row corresponds to the factor $\mathcal{S}(\mathbf{R})$, in this
case it is symmetric because our system consists of bosons. The advantage of
using $\Psi_{\rm S}(\mathbf{R})$ is that it is suitable for describing systems
with pairing. In particular, $\Psi_{\rm S}(\mathbf{R})$ describes well the
dimer-dimer problem.

Once the trial wave function has been chosen it is necessary to calculate its
gradient and Laplacian in order to implement the QMC method. The expressions
for the gradient and the Laplacian of $\Psi_{\rm T}(\mathbf{R})$ and
$\Psi_{\rm S}(\mathbf{R})$ can be found in the Appendices A and B,
respectively.

###  2.5 Quantum Monte Carlo Estimators

The aim of this section is to show how the ground-state properties are
computed in the Monte Carlo algorithm.

####  2.5.1 Pair Distribution Function

The pair distribution function $g(\mathbf{r}_{1},\mathbf{r}_{2})$ is
proportional to the probability of finding two particles at the positions
$\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, simultaneously. In coordinate
representation, for a homogeneous system $g(\mathbf{r}_{1},\mathbf{r}_{2})$ is
given by

$g(\mathbf{r}_{1},\mathbf{r}_{2})=\frac{N(N-1)}{n^{2}}\frac{\int|\Psi(\mathbf{R})|^{2}d\mathbf{r}_{3}\cdots
d\mathbf{r}_{N}}{\int|\Psi(\mathbf{R})|^{2}d\mathbf{R}},$ (2.77)

where $n$ is the density of the system. In a homogeneous system the pair
distribution function $g(\mathbf{r}_{1},\mathbf{r}_{2})$ depends only on the
relative distance $\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$, with this
assumption $g(\mathbf{r}_{1},\mathbf{r}_{2})$ becomes

$g(\mathbf{r})=\frac{N(N-1)}{n^{2}L^{d}}\frac{\int|\Psi(\mathbf{R})|^{2}\delta(\mathbf{r}_{12}-\mathbf{r})d\mathbf{R}}{\int|\Psi(\mathbf{R})|^{2}d\mathbf{R}},$
(2.78)

where $L$ is the size of the simulation box and $d$ is the dimensionality of
the system. To improve the efficiency of the calculation we sum over all pair
of particles

$g(\mathbf{r})=\frac{2}{nN}\frac{\int|\Psi(\mathbf{R})|^{2}\sum_{i<j}\delta(\mathbf{r}_{ij}-\mathbf{r})d\mathbf{R}}{\int|\Psi(\mathbf{R})|^{2}d\mathbf{R}},$
(2.79)

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$. In Monte Carlo, the
pair distribution function is determined by making a histogram of the distance
between all pair of particles in the system.

####  2.5.2 One-Body Density Matrix

For a homogeneous system described by the many-body wave function
$\Psi(\mathbf{r}_{1},\cdots,\mathbf{r}_{N})$ the one body density matrix
$\rho(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime})$ is defined as

$\rho(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime})=N\frac{\int
d\mathbf{r}_{2}\cdots
d\mathbf{r}_{N}\Psi^{*}(\mathbf{r}_{1},\mathbf{r}_{2},\cdots,\mathbf{r}_{N})\Psi(\mathbf{r}_{1}^{\prime},\mathbf{r}_{2},\cdots,\mathbf{r}_{N})}{\int
d\mathbf{r}_{1}\cdots
d\mathbf{r}_{N}|\Psi(\mathbf{r}_{1},\mathbf{r}_{2},\cdots,\mathbf{r}_{N})|^{2}}.$
(2.80)

In the VMC calculations, we sample the trial wave function $\Psi_{\rm
T}(\mathbf{R})$. Thus, a variational estimation of the one body density matrix
is given by

$\rho(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime})=N\frac{\int
d\mathbf{r}_{2}\cdots d\mathbf{r}_{N}\frac{\Psi_{\rm
T}^{*}(\mathbf{R})}{\Psi_{\rm T}^{*}(\mathbf{R}^{\prime})}|\Psi_{\rm
T}(\mathbf{R}^{\prime})|^{2}}{\int d\mathbf{R}|\Psi_{\rm
T}(\mathbf{R})|^{2}},$ (2.81)

where $\mathbf{R}=\\{\mathbf{r}_{1},\mathbf{r}_{2},\cdots,\mathbf{r}_{N}\\}$
and
$\mathbf{R}^{\prime}=\\{\mathbf{r}^{\prime}_{1},\mathbf{r}_{2},\cdots,\mathbf{r}_{N}\\}$.

Instead, in the DMC method we sample the mixed distribution
$f(\mathbf{R})=\Psi_{\rm T}(\mathbf{R})\Psi(\mathbf{R})$. Thus, a mixed
estimation of the one body density matrix is given by

$\rho(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime})=N\frac{\int
d\mathbf{r}_{2}\cdots d\mathbf{r}_{N}\frac{\Psi_{\rm
T}^{*}(\mathbf{R})}{\Psi_{\rm
T}^{*}(\mathbf{R}^{\prime})}f(\mathbf{R}^{\prime})}{\int
d\mathbf{R}f(\mathbf{R})}.$ (2.82)

For a homogeneous Bose system, the condensate fraction $N_{0}/N$ is obtained
from the asymptotic behavior of the one body density matrix

$\lim_{|\mathbf{r}-\mathbf{r^{\prime}}|\to\infty}\frac{\rho(\mathbf{r}_{1},\mathbf{r}_{1}^{\prime})}{n}=\frac{N_{0}}{N},$
(2.83)

where $N_{0}$ is the number of particles in the condensate.

####  2.5.3 Mixed Estimators and Extrapolation Technique

The expectation value of a given observable $\hat{O}$ is obtained from

$\langle\hat{O}\rangle=\frac{\langle\Psi|\hat{O}|\Psi\rangle}{\langle\Psi|\Psi\rangle},$
(2.84)

where $\Psi$ is the wave function of the system. In a VMC calculation, the
expectation values are evaluated with the trial wave function $\Psi_{\rm T}$,
therefore we obtain a variational estimator

$\langle\hat{O}\rangle_{\rm var}=\frac{\langle\Psi_{\rm T}|\hat{O}|\Psi_{\rm
T}\rangle}{\langle\Psi_{\rm T}|\Psi_{\rm T}\rangle}.$ (2.85)

A more precise estimator is obtained from a DMC calculation, where the
expectation values are sample for the mixed distribution $f=\Psi_{\rm T}\Psi$.
After long enough imaginary time propagation we have $f\approx\Psi_{\rm
T}\Phi_{0}$, and the expectation value is obtained from

$\langle\hat{O}\rangle_{\rm mix}=\frac{\langle\Psi_{\rm
T}|\hat{O}|\Phi_{0}\rangle}{\langle\Psi_{\rm T}|\Phi_{0}\rangle},$ (2.86)

where $\Phi_{0}$ is the ground-state wave function. The last equation is known
as the mixed estimator since it is calculated over two different states. The
Eq. (2.86) gives the exact expectation value for the Hamiltonian (i.e. the
calculation of the ground-state energy is exact) and for observables that
commute with it. In the case of operators that do not commute with $\hat{H}$,
the result obtained from Eq. (2.86) will be biased by $\Psi_{\rm T}$. In this
case, it is possible to improve the description by employing a first order
correction in $\Psi_{\rm T}$ using the extrapolation method. In this method,
one assumes that the difference between the trial wave function $\Psi_{\rm T}$
and the ground-state wave function $\Phi_{0}$ is small:
$\delta\Psi=\Phi_{0}-\Psi_{\rm T}$. The approximated value of the exact
estimator with a second-order error in $\delta\Psi$ can be written in two
forms

$\displaystyle\langle\hat{O}\rangle_{\rm ext_{1}}=$ $\displaystyle
2\langle\hat{O}\rangle_{\rm mix}-\langle\hat{O}\rangle_{\rm
var}+\mathcal{O}(\delta\Psi^{2}),$ (2.87)
$\displaystyle\langle\hat{O}\rangle_{\rm ext_{2}}=$
$\displaystyle\frac{\langle\hat{O}\rangle_{\rm
mix}^{2}}{\langle\hat{O}\rangle_{\rm var}}+\mathcal{O}(\delta\Psi^{2}).$

The main limitation of using the extrapolated estimators
$\langle\hat{O}\rangle_{\rm ext_{1}}$ and $\langle\hat{O}\rangle_{\rm
ext_{2}}$ is that they depend on the quality of trial wave function $\Psi_{\rm
T}$ used for importance sampling. However, it is usefull to have two different
estimators, as if they differ among themselves, the difference will show the
typical difference with the exact result.

####  2.5.4 Pure Estimators

To overcome the limitations of the extrapolation method, one can use forward
walking techniques or similar methods to calculate pure estimators for local
observables that do not commute with the Hamiltonian. The pure estimator of a
local observable $\hat{O}$ is given by

$\langle\hat{O}\rangle_{\rm
pure}=\frac{\langle\Phi_{0}|\hat{O}|\Phi_{0}\rangle}{\langle\Phi_{0}|\Phi_{0}\rangle}.$
(2.88)

The natural outcome of the DMC method is instead a mixed estimator Eq. (2.86),
which differs from Eq. (2.88) by the presence of the trial wave function
$\Psi_{\rm T}$ on one of the sides. Nevertheless, the pure estimator can be
related to the mixed one by reweighting the observable with the quotient
$\Phi_{0}/\Psi_{\rm T}$ calculated for the same coordinates as the local
observable

$\langle\hat{O}\rangle_{\rm pure}=\frac{\langle\Psi_{\rm
T}|\frac{\Phi_{0}}{\Psi_{\rm T}}\hat{O}|\Phi_{0}\rangle}{\langle\Psi_{\rm
T}|\frac{\Phi_{0}}{\Psi_{\rm
T}}|\Phi_{0}\rangle}=\left\langle\frac{\Phi_{0}}{\Psi_{\rm
T}}\hat{O}\right\rangle_{\rm mix}.$ (2.89)

According to Lie et. al. [54] the quotient $\Phi_{0}/\Psi_{\rm T}$ can be
computed from the asymptotic number of descendants of each of the walkers

$W(\mathbf{R})=\lim_{\tau\to\infty}n\left(\mathbf{R}(\tau)\right).$ (2.90)

Usually, in the Monte Carlo algorithm the local observables are calculated by
taking block averages. Each block consists of $M$ time steps or iterations.
Inside one of these blocks, after one iteration, when a walker is replicated,
we replicate its coordinates and its weight Eq. (2.90), and computed the
observables associated with it [55, 56]:

$\displaystyle O_{k,\alpha}$
$\displaystyle=\langle\hat{O}\left(\mathbf{R}_{k,\alpha}\right)\rangle_{\rm
mix},$ (2.91) $\displaystyle W_{k,\alpha}$
$\displaystyle=n\left(\mathbf{R}_{k,\alpha}\right),$

where $k$ is the time step index and $\alpha$ is an index over the number of
walkers $N_{w,k}$ After a block is completed, the estimator of the observable
is calculated as

$\langle\hat{O}\rangle^{\rm
block}=\frac{\sum_{k=0}^{M}\sum_{\alpha=1}^{N_{w,k}}W_{k,\alpha}O_{k,\alpha}}{\sum_{k=0}^{M}\sum_{\alpha=1}^{N_{w,k}}W_{k,\alpha}}.$
(2.92)

After $N_{block}$ blocks, the pure estimator is given by

$\langle\hat{O}\rangle_{\rm pure}=\frac{1}{N_{\rm block}}\sum_{j=1}^{N_{\rm
block}}\langle\hat{O}\rangle^{j}.$ (2.93)

The pure estimator depends on the size $M$ of a block. $M$ has to be large
enough to reach the asymptotic regime given by Eq. (2.90).

Although the calculation of the pure estimators for the potential energy,
density profile, pair distribution function, static structure factor, and
other correlation functions are routinely done to our best knowledge, the
calculation of pure coordinates never has been done. We found it convenient to
store the coordinates of the walker and replicate them during the branching
process. After long enough propagation time, the pure coordinates are stored
and at the end of the simulation, we have a large number of pure coordinates.
At this point an average of a local observable over them automatically becomes
pure. In particular, we find this trick to be very flexible and especially
useful for the pure estimation of all sorts of complicated correlation
functions common for few-body analysis and often involving calculations of
hyperradius. For example, Fig. 5.3, Fig. 5.4 and Fig. 5.5, were obtained by
this method.

##  Chapter 3 One-Dimensional Three-Boson Problem with Two- and Three-Body
Interactions

In this chapter, we study the three-boson problem with contact two- and three-
body interactions in one dimension. By using the diffusion Monte Carlo
technique we calculate the binding energy of two and three dimers formed in a
Bose-Bose or Fermi-Bose mixture with attractive interspecies and repulsive
intraspecies interactions. Combining these results with a three-body theory
[57], we extract the three-dimer scattering length close to the dimer-dimer
zero crossing. In both considered cases the three-dimer interaction turns out
to be repulsive. Our results constitute a concrete proposal for obtaining a
one-dimensional gas with a pure three-body repulsion.

###  3.1 Introduction

The one-dimensional $N$-boson problem with the two-body contact interaction
$g_{2}\delta(x)$ is exactly solvable. Lieb and Liniger [58] have shown that
for $g_{2}>0$ the system is in the gas phase with positive compressibility.
McGuire [59] has demonstrated that for $g_{2}<0$ the ground state is a soliton
with the chemical potential diverging with $N$. In the case $N=\infty$ the
limits $g_{2}\rightarrow+0$ and $g_{2}\rightarrow-0$ are manifestly different:
The former corresponds to an ideal gas whereas the latter corresponds to
collapse. Accordingly, the behavior of a realistic one- or quasi-one-
dimensional system close to the two-body zero crossing strongly depends on
higher-order terms not included in the Lieb-Liniger or McGuire zero-range
models. Sekino and Nishida [60] have considered one-dimensional bosons with a
pure zero-range three-body attraction and found that the ground state of the
system is a droplet with the binding energy exponentially increasing with $N$,
which also means collapse in the thermodynamic limit. In Ref. [61], the
authors have argued that in a sufficiently dilute regime the three-body
interaction is effectively repulsive, providing a mechanical stabilization
against collapse for $g_{2}<0$. The competition between the two-body
attraction and three-body repulsion leads to a dilute liquid state similar to
the one discussed by Bulgac [62] in three dimensions.

The three-body scattering in one dimension is kinematically equivalent to a
two-dimensional two-body scattering [63, 60]. Therefore, the corresponding
interaction shift depends logarithmically on the product of the scattering
momentum and three-body scattering length $a_{3}$. An important consequence of
this fact is that, in contrast to higher dimensions, the one-dimensional
three-body interaction can become noticeable even if $a_{3}$ is exponentially
small compared to the mean interparticle distance. Therefore, three-body
effects can be studied in the universal dilute regime essentially in any one-
dimensional system that preserves a finite residual three-body interaction
close to a two-body zero crossing. Universality means that the effective-range
effects are exponentially small and the relevant interaction parameters are
the two- and three-body scattering lengths $a_{2}$ and $a_{3}$, respectively.

In this chapter, we consider a two-component Bose-Bose mixture with attractive
interspecies and repulsive intraspecies interactions. In this system, the
interspecies attraction binds atoms into dimers while the dimer-dimer
interaction is tunable by changing the intraspecies repulsion [61]. Analytical
predictions [57] are complemented by diffusion Monte Carlo calculations of the
hexamer energy, permitting to determine the three-dimer scattering length
close to the dimer-dimer zero crossing. We perform this procedure for equal
intraspecies coupling constants and in the case where their ratio is infinite.
In the latter limit, one of the components is in the Tonks-Girardeau regime
and the system is equivalent to a Fermi-Bose mixture. We find that the three-
dimer interaction is repulsive in both cases.

###  3.2 The System

In Ref. [57], the authors considered a one-dimensional system of three bosons
of mass $m$ interacting via contact two- and three-body forces characterized
by the scattering lengths $a_{2}$ and $a_{3}$, respectively. They obtained the
following analytical expression for the ground and excited trimer energies

$\ln\frac{a_{3}\kappa
e^{\gamma}}{a_{2}}=\frac{2}{\kappa^{2}-1}\left[\frac{\pi}{3\sqrt{3}}+\frac{3\kappa^{2}-1}{\sqrt{4\kappa^{2}-1}}\arctan\sqrt{\frac{2\kappa+1}{2\kappa-1}}\right],$
(3.1)

where $\kappa=\sqrt{-mE}a_{2}/(2\hbar)$ and $\gamma=0.577$ is Euler’s
constant. The Eq. (3.1) relates the trimer binding energy $E=E_{3}<0$ with
$a_{2}$ and $a_{3}$. Considering the dimer binding energy as
$|E_{2}|=\hbar^{2}/ma_{2}^{2}$ we obtain the following relation
$E_{3}/E_{2}=4\kappa^{2}$. In the following we are going to use Eq. (3.1) to
extract the three-dimer scattering length close to the dimer-dimer zero
crossing.

Systems where two- and three-body effective interactions can be controlled
independently are difficult to produce or engineer (see [63] and references
therein). We now discuss a model tunable to the regime of pure three-body
repulsion. Namely, we consider a mixture of one-dimensional pointlike bosons A
and B of unit mass characterized by the coupling constants

$g_{{\text{A}}{\text{B}}}=-\frac{2\hbar^{2}}{ma_{{\text{A}}{\text{B}}}}<0,$
(3.2)

for the interspecies attraction and

$g_{\sigma\sigma}=-\frac{2\hbar^{2}}{ma_{\sigma\sigma}}>0,$ (3.3)

for the intraspecies repulsions. The interspecies attraction leads to the
formation of AB dimers of size $a_{{\text{A}}{\text{B}}}$ and energy

$E_{{\text{A}}{\text{B}}}=-\frac{\hbar^{2}}{ma_{{\text{A}}{\text{B}}}^{2}}.$
(3.4)

One can show [61] that the two-dimer interaction changes from attractive to
repulsive with increasing $g_{\sigma\sigma}$. In particular, the two-dimer
zero crossing is predicted to take place for the Bose-Bose (BB) case at

$g_{{\text{A}}{\text{A}}}=g_{{\text{B}}{\text{B}}}=2.2|g_{{\text{A}}{\text{B}}}|,$
(3.5)

and fo the Fermi-Bose (FB) case at

$g_{{\text{A}}{\text{A}}}=0.575|g_{{\text{A}}{\text{B}}}|,$ (3.6)

if $g_{{\text{B}}{\text{B}}}=\infty$.

Here, we consider three such dimers and characterize their three-dimer
interaction by calculating the hexamer energy
$E_{{\text{A}}{\text{A}}{\text{A}}{\text{B}}{\text{B}}{\text{B}}}$ and by
comparing it with the tetramer energy
$E_{{\text{A}}{\text{A}}{\text{B}}{\text{B}}}$ on the attractive side of the
two-dimer zero crossing where the tetramer exists. The idea is that
sufficiently close to this crossing the dimers behave as pointlike particles
weakly bound to each other. One can then extract the three-dimer scattering
length $a_{3}$ from the zero-range three-boson formalism [Eq. (3.1)], with
$m\to 2m$

$\displaystyle E_{2}$
$\displaystyle=E_{{\text{A}}{\text{A}}{\text{B}}{\text{B}}}-2E_{{\text{A}}{\text{B}}},$
(3.7) $\displaystyle E_{3}$
$\displaystyle=E_{{\text{A}}{\text{A}}{\text{A}}{\text{B}}{\text{B}}{\text{B}}}-3E_{{\text{A}}{\text{B}}},$

and using the asymptotic expression for the dimer-dimer scattering length

$E_{2}=-\frac{\hbar^{2}}{2ma_{2}^{2}},$ (3.8)

we obtain

$a_{2}=\frac{\hbar}{\sqrt{2m|E_{2}|}}.$ (3.9)

###  3.3 Details of the Methods

In order to calculate $E_{2}$ and $E_{3}$, we resort to the diffusion Monte
Carlo (DMC) technique,, which was explained in Chapter 2. The importance
sampling is used to reduce the statistical noise and also to impose the Bethe-
Peierls boundary conditions stemming from the $\delta$-function interactions.
We construct the guiding wave function $\psi_{\rm T}$ in the pair-product form

$\psi_{\rm
T}=\prod\limits_{i<j}f^{{\text{A}}{\text{A}}}(x_{ij}^{{\text{A}}{\text{A}}})\prod\limits_{i<j}f^{{\text{B}}{\text{B}}}(x_{ij}^{{\text{B}}{\text{B}}})\prod\limits_{i,j}f^{{\text{A}}{\text{B}}}(x_{ij}^{{\text{A}}{\text{B}}})\;,$
(3.10)

where $x_{ij}^{\sigma\sigma^{\prime}}=x_{i}^{\sigma}-x_{j}^{\sigma^{\prime}}$
is the distance between particles $i$ and $j$ of components $\sigma$ and
$\sigma^{\prime}$, respectively. The intercomponent correlations are governed
by the dimer wave function

$f^{{\text{A}}{\text{B}}}(x)=\exp\left(-\frac{|x|}{a_{{\text{A}}{\text{B}}}}\right),$
(3.11)

and the intracomponent terms are

$f^{\sigma\sigma}(x)=\sinh\left(\frac{|x|}{a_{{\text{A}}{\text{B}}}}-\frac{|x|}{2a_{\rm
dd}}\right)-\left(\frac{a_{\sigma\sigma}}{a_{{\text{A}}{\text{B}}}}-\frac{a_{\sigma\sigma}}{2a_{\rm
dd}}\right).$ (3.12)

These functions satisfy the Bethe-Peierls boundary conditions,

$\frac{\partial f^{\sigma\sigma^{\prime}}(x)}{\partial
x}|_{x=+0}=-\frac{f^{\sigma\sigma^{\prime}}(0)}{a_{\sigma\sigma^{\prime}}},$
(3.13)

which, because of the product form, also ensures the correct behavior of the
total guiding function $\psi_{\rm T}$ at any two-body coincidence. At the same
time, the long-distance behavior of $f^{\sigma\sigma}(x)$ is chosen such that
$\psi_{\rm T}$ allows dimers to be at distances larger than their size. When
the distance $x$ between pairs $\\{x_{1}^{\text{A}},x_{1}^{\text{B}}\\}$ and
$\\{x_{2}^{\text{A}},x_{2}^{\text{B}}\\}$ is much larger than the dimer size
$a_{{\text{A}}{\text{B}}}$, Eq. (3.10) reduces to

$\psi_{\rm T}\propto
f^{{\text{A}}{\text{B}}}(x_{11}^{{\text{A}}{\text{B}}})f^{{\text{A}}{\text{B}}}(x_{22}^{{\text{A}}{\text{B}}})\exp\left(-\frac{|x|}{a_{\rm
dd}}\right).$ (3.14)

For $a_{\rm dd}\gg a_{{\text{A}}{\text{B}}}$, this wave function describes two
dimers weakly-bound to each other. It can be noted that the choice of the same
spin Jastrow terms might seem rather unusual as $f^{\sigma\sigma}$ become
exponentially large at large distances due to divergence of $\sinh(x)$
function. This divergence is cured by multiplication of exponentially small
opposite-spin Jatsrow terms $f^{{\text{A}}{\text{B}}}(x)$. Thus, the use of
the $\sinh(x)$ function allows both to impose the physically correct long-
range properties, Eq. (3.14) and the correct Bethe-Peierls boundary condition
at short distances where the expansion $\sinh(x)\to x$ results in
$f^{\sigma\sigma}\propto|x|-a^{\sigma\sigma}$.

While $a_{\sigma\sigma^{\prime}}$ are fixed by the Hamiltonian, we treat
$a_{\rm dd}$ as a free parameter in Eq. (3.10). Close to the dimer-dimer zero
crossing $a_{\rm dd}\approx a_{2}$ and this parameter is related self-
consistently to the tetramer energy while far from the crossing its value is
optimized according to the variational principle. It is useful to mention that
in case FB, where $a_{{\text{B}}{\text{B}}}=0$, the B component is in the
Tonks-Girardeau limit and can be mapped to ideal fermions by Girardeau’s
mapping [64]. Replacing $|x|$ by $x$ in the definition of
$f^{{\text{B}}{\text{B}}}(x)$ makes $\psi_{\rm T}$ antisymmetric with respect
to permutations of B coordinates.

###  3.4 Results

####  3.4.1 Binding Energies

In Fig. 3.1 (left) we show the tetramer energy
$E_{{\text{A}}{\text{A}}{\text{B}}{\text{B}}}$ as a function of the ratio
$g_{{\text{A}}{\text{A}}}/|g_{{\text{A}}{\text{B}}}|$ for the Bose-Fermi and
Bose-Bose mixtures. The thresholds for binding are shown by arrows.

In Fig. 3.1 (right) we show the hexamer energy
$E_{{\text{A}}{\text{A}}{\text{A}}{\text{B}}{\text{B}}{\text{B}}}$ as a
function of the ratio $g_{{\text{A}}{\text{A}}}/|g_{{\text{A}}{\text{B}}}|$
for the Bose-Fermi and Bose-Bose mixtures. We find that for sufficiently
strong intercomponent repulsion (larger
$g_{{\text{A}}{\text{A}}}/|g_{{\text{A}}{\text{B}}}|$) the hexamer gets
unbound, first for the Bose-Fermi case and then for the Bose-Bose mixture.

Figure 3.1: Tetramer $E_{{\text{A}}{\text{A}}{\text{B}}{\text{B}}}$ (left) and
hexamer $E_{{\text{A}}{\text{A}}{\text{A}}{\text{B}}{\text{B}}{\text{B}}}$
(right) energies in units of the dimer energy $|E_{{\text{A}}{\text{B}}}|$ for
Bose-Fermi and Bose-Bose mixtures as a function of the ratio
$g_{{\text{A}}{\text{A}}}/|g_{{\text{A}}{\text{B}}}|$. The arrows show the
positions of the thresholds for binding.

####  3.4.2 Threshold Determination

In Fig. 3.2 we show the numerical threshold determination for the tetramer and
hexamer for the Bose-Fermi and Bose-Bose mixtures. Our numerical results for
the tetramer threshold values are consistent with the predictions of Ref. [61]
(Eq. 3.5 and Eq. 3.6). We find that the hexamer threshold for the Bose-Fermi
mixture is located at
$g_{{\text{A}}{\text{A}}}/|g_{{\text{A}}{\text{B}}}|\approx 0.575$. In the
case of the Bose-Bose mixture, the hexamer threshold occurs at
$g_{{\text{A}}{\text{A}}}/|g_{{\text{A}}{\text{B}}}|\approx 2.2$. In both
cases, the hexamer thresholds coincide with the tetramer thresholds.

Figure 3.2: Fitting procedure for the determination of the threshold values of
the AABB tetramer and AABBB hexamer for Bose-Fermi and Bose-Bose mixtures. The
dashed curves correspond to the fit
$c_{1}\left(\frac{g_{\mathrm{AA}}}{|g_{\mathrm{AB}}|}-\frac{g_{\mathrm{AA}}^{c}}{|g_{\mathrm{AB}}^{c}|}\right)+c_{2}\left(\frac{g_{\mathrm{AA}}}{|g_{\mathrm{AB}}|}-\frac{g_{\mathrm{AA}}^{c}}{|g_{\mathrm{AB}}^{c}|}\right)^{2}$.

####  3.4.3 Three-Dimer Repulsion

In Fig. 3.3, we show $E_{3}/|E_{2}|$ for cases BB (red squares) FB (blue
circles) as a function of $\delta=1/\ln(\sqrt{2m|E_{2}|}a_{3}/\hbar)$ along
with the prediction of Eq. (3.1) (solid black). The quantity $a_{3}$ is a
fitting parameter to the DMC results; changing it essentially shifts the data
horizontally. We clearly see that in both cases the three-dimer interaction is
repulsive since $E_{3}/|E_{2}|$ is above the McGuire trimer limit [59] (dash-
dotted line). For rightmost data points the hexamer is about ten times larger
than the dimer and the data align with the universal zero-range analytics. For
the other points we observe significant effective range effects related to the
finite size of the dimer. In the universal limit $a_{{\text{A}}{\text{B}}}\ll
a_{2}$, the leading effective-range correction to the ratio $E_{3}/|E_{2}|$ is
expected to be proportional to $a_{{\text{A}}{\text{B}}}/a_{2}\propto
e^{1/\delta}$ [61]. Indeed, adding the term $Ce^{1/\delta}$ to the zero-range
prediction well explains deviations of our results from the universal curve
and we have checked that other exponents do not work that well. We thus treat
$a_{3}$ and $C$ as fitting parameters; in case BB we obtain
$a_{3}=0.01a_{{\text{A}}{\text{B}}}$ and in case FB
$a_{3}=0.03a_{{\text{A}}{\text{B}}}$. Both cases are fit with $C=-100$ (dashed
curve in Fig. 3.3). We emphasize that we are dealing with the true ground
state of three dimers. The lower “attractive” state formally existing for
these values of $a_{2}$ and $a_{3}$ in the zero-range model is an artifact
since it does not satisfy the zero-range applicability condition. The three-
dimer interaction is an effective finite-range repulsion which supports no
bound states.

Figure 3.3:  $E_{3}/|E_{2}|$ vs $1/\ln(\sqrt{2m|E_{2}|}a_{3}/\hbar)$ for one-
dimensional dimers. Here $E_{2}$ and $E_{3}$ are the tetramer and hexamer
energies measured relative to the two- and three-dimer thresholds,
respectively. The solid curve is the prediction of Eq. (3.1) and the dashed
curve is a fit, which includes finite-dimer-size effects into account. The
dash-dotted line is the McGuire result $E_{3}=4E_{2}$ for three pointlike
bosons with no three-body interaction. The red squares are the DMC data for
case BB plotted using $a_{3}=0.01a_{{\text{A}}{\text{B}}}$ and the blue
circles stand for case FB with $a_{3}=0.03a_{{\text{A}}{\text{B}}}$. The error
bars are larger in the latter case because of the larger statistical noise
induced by the nodal surface imposed by the Fermi statistics.

###  3.5 Summary

In conclusion, we argue that since in one dimension the three-body energy
correction scales logarithmically with the three-body scattering length
$a_{3}$, three-body effects are observable even for exponentially small
$a_{3}$, which significantly simplifies the task of engineering three-body-
interacting systems in one dimension. We demonstrate that Bose-Bose or Fermi-
Bose dimers, previously shown to be tunable to the dimer-dimer zero crossing,
exhibit a noticeable three-dimer repulsion. We can now be certain that the
ground state of many such dimers slightly below the dimer-dimer zero crossing
is a liquid in which the two-body attraction is compensated by the three-body
repulsion [62, 61].

Our results have implications for quasi-one-dimensional mixtures. We mention
particularly the 40K-41K Fermi-Bose mixture which emerges as a suitable
candidate for exploring the liquid state of fermionic dimers. Here the
intraspecies 41K-41K background interaction is weakly repulsive (the triplet
41K-41K scattering length equals 3.2nm [65]) and the interspecies one features
a wide Feshbach resonance at 540G [66]. Let us identify A with 40K, B with
41K, and assume the radial oscillator length $l_{0}=56$nm, which corresponds
to the confinement frequency $2\pi\times 80$kHz. Under these conditions the
effective coupling constants equal $g_{\sigma\sigma^{\prime}}\approx
2a_{\sigma\sigma^{\prime}}^{\rm(3D)}/l_{0}^{2}$ [67] and the dimer-dimer zero
crossing at $g_{{\text{A}}{\text{A}}}=0.575|g_{{\text{A}}{\text{B}}}|$ is
realized for the three-dimensional scattering lengths
$a_{{\text{B}}{\text{B}}}^{\rm(3D)}\approx 3.2$nm and
$a_{{\text{A}}{\text{B}}}^{\rm(3D)}\approx-5.6$nm. The dimer size is then
$\approx 560$nm and dimer binding energy corresponds to $\approx 2\pi\times
800$Hz placing the system in the one-dimensional regime. For the rightmost
(next to rightmost) blue circle in Fig. 3.3, the tetramer is approximately 20
(10) times larger than the dimer and 800 (200) times less bound. Moving left
in this figure is realized by increasing $|a_{{\text{A}}{\text{B}}}^{\rm 3D}|$
and thus getting deeper in the region
$g_{{\text{A}}{\text{A}}}<0.575|g_{{\text{A}}{\text{B}}}|$. Note, however,
that this also pushes the system out of the one-dimensional regime and effects
of transversal modes [68, 69, 70] become important.

##  Chapter 4 Few-Body Bound States of Two-Dimensional Bosons

In this chapter, we study clusters of the type ANBM with $N\leq M\leq 3$ in a
two-dimensional mixture of A and B bosons, with attractive AB and equally
repulsive AA and BB interactions. In order to check universal aspects of the
problem, we choose two very different models: dipolar bosons in a bilayer
geometry (this work) and particles interacting via separable Gaussian
potentials (reported in Ref. [71]). We find that all the considered clusters
are bound and that their energies are universal functions of the scattering
lengths $a_{{\text{A}}{\text{B}}}$ and
$a_{{\text{A}}{\text{A}}}=a_{{\text{B}}{\text{B}}}$, for sufficiently large
attraction-to-repulsion ratios
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$. When
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$ decreases below $\approx
10$, the dimer-dimer interaction changes from attractive to repulsive and the
population-balanced AABB and AAABBB clusters break into AB dimers. Calculating
the AAABBB hexamer energy just below this threshold, we find an effective
three-dimer repulsion which may have important implications for the many-body
problem, particularly for observing liquid and supersolid states of dipolar
dimers in the bilayer geometry. The population-imbalanced ABB trimer, ABBB
tetramer, and AABBB pentamer remain bound beyond the dimer-dimer threshold. In
the dipolar model, they break up at $a_{{\text{A}}{\text{B}}}\approx
2a_{{\text{B}}{\text{B}}}$ where the atom-dimer interaction switches to
repulsion. The work presented in this chapter was a collaboration [71]. I did
the calculations for the dipolar clusters.

###  4.1 Introduction

Recent experiments on dilute quantum droplets in dipolar bosonic gases [9, 10,
11, 12] and in Bose-Bose mixtures [6, 7, 8] with competing interactions have
exposed the important role of beyond-mean-field effects in weakly-interacting
systems. A natural strategy to boost these effects and enhance exotic
behaviors is to make the interactions stronger while keeping the attraction-
repulsion balance for mechanical stability. The most straightforward way of
getting into this regime is to increase the gas parameter $na_{s}^{3}$.
However, this leads to enhanced three-body losses which results in very short
lifetimes (as it has been observed in experiments [9, 10, 11, 12, 6, 7, 8]).
Nevertheless, this regime is achievable in reduced geometries. It has been
shown that a one-dimensional Bose-Bose mixture with strongly-attractive
interspecies interaction becomes dimerized and, by increasing the intraspecies
repulsion, the dimer-dimer interaction can be tuned from attractive to
repulsive [72]. Then, an effective three-dimer repulsion has been found in
this system and predicted to stabilize a liquid phase of attractive dimers
[57].

In two dimensions, a particularly interesting realization of such a strongly-
interacting, tunable, and long-lived Bose-Bose mixture is a system of dipolar
bosons confined to a bilayer geometry [73, 74, 75]. When the dipoles are
oriented perpendicularly to the plane, there is a competing effect between
repulsive intralayer and partially attractive interlayer interactions,
interesting from the viewpoint of liquid formation. In addition, the quasi-
long range character of the dipolar interaction can produce the rotonization
of its spectrum and a supersolid behavior [76, 77, 78, 79, 80, 81, 82, 83,
84], formation of a crystal phase [30, 31], and a pair superfluid [32, 33, 34]
(see also lattice calculations of Ref. [85]). A peculiar feature of bilayer
model is the vanishing Born integral for the interlayer interaction [86],

$\int V_{{\text{A}}{\text{B}}}(\rho)d^{2}\rho=0,$ (4.1)

which has led to controversial claims about the existence of a two-body bound
state [87] till it has finally been established that this bound state always
exists, although its energy can be exponentially small [88, 89, 90, 91, 92],
consistently with Ref. [93]. Interestingly, a similar controversy seems to
continue at the few-body level; it has been claimed [94] that the repulsive
dipolar tails will never allow for three- or four-body bound states in this
geometry.

In this chapter, we investigate few-body bound states in a two-dimensional
mass-balanced mixture of A and B bosons with two types of interactions
characterized by the two-dimensional scattering lengths $a_{\rm AB}$ and
$a_{\rm AA}=a_{\rm BB}$. The first case corresponds to the bilayer of dipoles
discussed above and, in the second, the interactions were modeled by non-local
(separable) finite-range Gaussian potentials [71]. By using the diffusion
Monte Carlo (DMC) technique in the first case, and the Stochastic Variational
Method (SVM) in the second, we find that for sufficiently weak BB repulsion
compared to the AB attraction, $a_{\rm AB}\gg a_{\rm BB}$, all clusters of the
type ANBM with $1\leq N\leq M\leq 3$ are bound. We then locate thresholds for
their unbinding with decreasing
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$. By looking at the AAABBB
hexamer energy close to the corresponding threshold, we discover an effective
three-dimer repulsion, which can stabilize interesting many-body phases.

###  4.2 The Hamiltonian

The Hamiltonian of the system is

$\displaystyle\hat{H}=$
$\displaystyle-\frac{\hbar^{2}}{2m}\sum_{i=1}^{N}\nabla^{2}_{i}-\frac{\hbar^{2}}{2m}\sum_{\alpha=1}^{M}\nabla_{\alpha}^{2}+\sum_{i<j}\hat{V}_{{\text{A}}{\text{A}}}(r_{ij})$
(4.2)
$\displaystyle+\sum_{\alpha<\beta}\hat{V}_{{\text{B}}{\text{B}}}(r_{\alpha\beta})+\sum_{i,\alpha}\hat{V}_{{\text{A}}{\text{B}}}(r_{i\alpha})\
,$

where the two-dimensional vectors ${\bf r}_{i}$ and ${\bf r}_{\alpha}$ denote
particle positions of species A and B containing, respectively, $N$ and $M$
atoms, $\hat{V}_{{\text{A}}{\text{B}}}$ and
$\hat{V}_{{\text{A}}{\text{A}}}=\hat{V}_{{\text{B}}{\text{B}}}$ are the
interspecies and intraspecies interaction potentials, and $m$ is the mass of
each particle. For the bilayer setup, we have

$V_{{\text{A}}{\text{A}}}(r)=V_{{\text{B}}{\text{B}}}(r)=\frac{d^{2}}{r^{3}},$
(4.3)

and

$V_{{\text{A}}{\text{B}}}(r)=\frac{d^{2}(r^{2}-2h^{2})}{(r^{2}+h^{2})^{5/2}},$
(4.4)

where $d$ is the dipole moment and $h$ is the distance between the layers.
Dipoles are aligned perpendicularly to the layers and there is no interlayer
tunneling. The potential $V_{{\text{B}}{\text{B}}}(r)$ is purely repulsive and
is characterized by the $h$-independent scattering length
$a_{{\text{B}}{\text{B}}}=e^{2\gamma}r_{0}$ [95], where $\gamma\approx 0.577$
is the Euler constant and $r_{0}=md^{2}/\hbar^{2}$ is the dipolar length. The
interlayer potential $V_{{\text{A}}{\text{B}}}(r)$ always supports at least
one dimer state. Its energy reported in the Fig. 4.1 diverges for $h\to 0$ and
exponentially vanishes in the opposite limit [88, 89, 90, 91]. The scattering
length $a_{{\text{A}}{\text{B}}}$, which is a function of $r_{0}$ and $h$, is
$\sim a_{{\text{B}}{\text{B}}}\sim r_{0}$ for $h\sim r_{0}$, and exponentially
large for $h\gg r_{0}$. In the following, we parametrize the system by
specifying $a_{{\text{B}}{\text{B}}}$ and $a_{{\text{A}}{\text{B}}}$ rather
than $h$ and $r_{0}$.

Figure 4.1: Dimer energy $E_{{\text{A}}{\text{B}}}$ in units of
$\hbar^{2}/ma_{{\text{B}}{\text{B}}}^{2}$ for Gaussian (curve) and dipolar
(symbols) potentials as a function of the scattering length ratio
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$.

In the more academic case of Gaussian interactions, the following potential
was used [71]

$\hat{V}_{{\text{A}}{\text{B}}}(r_{i\alpha})\psi(r_{i\alpha})=\int
V_{{\text{A}}{\text{B}}}(r_{i\alpha},r^{\prime}_{i\alpha})\psi(r^{\prime}_{i\alpha})d^{2}r^{\prime}_{i\alpha},$
(4.5)

and similarly for $V_{{\text{A}}{\text{A}}}$ and $V_{{\text{B}}{\text{B}}}$,
where

$\displaystyle V_{\sigma\sigma^{\prime}}(r,r^{\prime})$
$\displaystyle=C_{\sigma\sigma^{\prime}}G_{\xi}(r)G_{\xi}(r^{\prime}),$ (4.6)
$\displaystyle G_{\xi}(r)$
$\displaystyle=(2\pi\xi^{2})^{-1}\exp(-r^{2}/2\xi^{2}),$

and $\xi$ is the characteristic range of the potential. An advantage of this
non-local potential is that the two-body problem can be solved analytically,
giving

$C_{\sigma\sigma^{\prime}}^{-1}=\frac{m}{4\pi\hbar^{2}}\left[2\ln\frac{2\xi}{a_{\sigma\sigma^{\prime}}}-\gamma\right].$
(4.7)

In the following, the ratio is varied
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$, with
$a_{{\text{B}}{\text{B}}}=1.4\xi$ fixed. Note that the available ratio is
limited to $a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}>1.1$.

###  4.3 Details of the Methods

In order to calculate the energies of the different few-body clusters with
dipolar interactions we use the diffusion Monte Carlo (DMC) method (see
Chapter 2), which leads to the exact ground-state energy of the system, within
a statistical error. This stochastic technique solves the Schrödinger equation
in imaginary time using a trial wave function for importance sampling. We
choose it to be

$\displaystyle\Psi_{\rm S}(\mathbf{r}_{1},\dots,\mathbf{r}_{N+M})$
$\displaystyle=\prod_{i<j}^{N}f_{{\text{A}}{\text{A}}}(r_{ij})\prod_{\alpha<\beta}^{M}f_{{\text{B}}{\text{B}}}(r_{\alpha\beta})$
(4.8)
$\displaystyle\times\left[\prod_{i=1}^{N}\sum_{\alpha=1}^{M}f_{{\text{A}}{\text{B}}}(r_{i\alpha})+\prod_{\alpha=1}^{M}\sum_{i=1}^{N}f_{{\text{A}}{\text{B}}}(r_{i\alpha})\right],$

which takes into account a possible formation of AB dimers.

Figure 4.2: Intraspecies $f_{{\text{A}}{\text{A}}}(r)$ and interspecies
$f_{{\text{A}}{\text{B}}}(r)$ wave functions.

The intraspecies Jastrow factors are chosen as the zero-energy two-body
scattering solution,

$f_{{\text{A}}{\text{A}}}(r)=f_{{\text{B}}{\text{B}}}(r)=K_{0}(2\sqrt{r_{0}/r}),$
(4.9)

with $K_{0}$ the modified Bessel function. The interspecies interactions are
described by the dimer wave function $f_{{\text{A}}{\text{B}}}(r)$ up to
$R_{0}$, calculated numerically. The variational parameter $R_{0}$ is chosen
to be large enough that for distances larger than $R_{0}$ we neglected the
dipolar potential and took the free scattering solution
$f_{{\text{A}}{\text{B}}}(r)=CK_{0}(\sqrt{-mE_{{\text{A}}{\text{B}}}}r/\hbar)$.
We impose continuity of the logarithmic derivative at the matching point
$R_{0}$, this condition yields to the following equality

$\frac{f^{{}^{\prime}}_{{\text{A}}{\text{B}}}(R_{0})}{f_{{\text{A}}{\text{B}}}(R_{0})}=-\frac{\sqrt{-mE_{{\text{A}}{\text{B}}}}}{\hbar}\frac{K_{1}(\sqrt{-mE_{{\text{A}}{\text{B}}}}R_{0}/\hbar)}{K_{0}(\sqrt{-mE_{{\text{A}}{\text{B}}}}R_{0}/\hbar)}.$
(4.10)

In Fig. 4.2 we show the intraspecies $f_{{\text{A}}{\text{A}}}(r)$ and
interspecies $f_{{\text{A}}{\text{B}}}(r)$ wave functions.

In the Gaussian model, the stochastic variational method (SVM) was used.
Details of the SVM method can be found in Refs. [96, 97].

###  4.4 Results

####  4.4.1 Binding Energies

We first discuss the limit of very large $a_{{\text{A}}{\text{B}}}$ (large
dimer size) when the interaction range and the intraspecies interactions can
be neglected. In this case, the problem can be treated in the zero-range
approximation giving for the ABB trimer $E_{\rm
ABB}^{a_{{\text{B}}{\text{B}}}=0}=2.39E_{\rm AB}$ [98, 99, 100] and for the
tetramers $E_{\rm ABBB}^{a_{{\text{B}}{\text{B}}}=0}=4.1E_{\rm AB}$ and
$E_{\rm AABB}^{a_{{\text{B}}{\text{B}}}=0}=10.6E_{\rm AB}$ [98]. The other
ANBM clusters (with $1\leq N\leq M\leq 3$) are also bound in absence of the
intraspecies repulsion. In Ref. [101], the authors calculated their binding
energies (and they also updated the energies of smaller clusters), which are
reported in Table 4.1.

ANBM | $E^{a_{{\text{B}}{\text{B}}}=0}_{{\text{A}}_{N}{\text{B}}_{M}}/E_{{\text{A}}{\text{B}}}$
---|---
ABB | 2.3896(1)
ABBB | 4.1364(2)
AABB | 10.690(2)
AABBB | 28.282(5)
AAABBB | 104.01(5)
Table 4.1: Energies of ${\rm A}_{N}{\rm B}_{M}$ clusters in units of the dimer
energy $E_{{\text{A}}{\text{B}}}$ in absence of the intraspecies repulsion
$a_{{\text{B}}{\text{B}}}=0$ [101]. The number between parenthesis indicates
the error in the last digits of the corresponding value.

The intraspecies repulsion shifts the cluster energies upwards as has been
seen for the ABB trimer [102, 103] and for the ABBB tetramer [103]. In Fig.
4.3, we report the energies of these and bigger clusters for both the dipolar
and Gaussian interactions. Note that, even for the weakest BB repulsion shown
in this figure ($a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}=200$), the
clusters are significantly less bound compared to the case of no BB repulsion.
This happens since the small parameter that controls the weakness of the
intraspecies interaction relative to the interspecies one is actually
$\lambda=1/\ln(a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}})\ll 1$. By
contrast, effective-range corrections contain powers of $r_{0}\sqrt{mE}/\hbar$
or $\xi\sqrt{mE}/\hbar$ for dipolar or Gaussian interactions, respectively,
which are exponentially small in terms of $\lambda$. This explains why the two
interaction models lead to almost indistinguishable results for large
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$.

Figure 4.3: Energies of ${\rm A}_{N}{\rm B}_{M}$ clusters in units of the
dimer energy $E_{{\text{A}}{\text{B}}}$ for Gaussian (curves) and dipolar
(symbols) potentials as a function of the scattering length ratio
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$. The arrows show the
positions of the thresholds for binding in the dipolar case. Figure 4.4:
Binding energies of the few-body clusters $E_{{\rm A}_{N}{\rm B}_{M}}-NE_{\rm
AB}$, in units of $\hbar^{2}/ma_{{\text{B}}{\text{B}}}^{2}$, versus
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$, for Gaussian (curves) and
dipolar (symbols) potentials. The arrows show the positions of the thresholds
for binding in the dipolar case. The uppper axis is labeled in dipolar units.

We find that for sufficiently strong intraspecies repulsion (smaller
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$) the trimer and all higher
clusters get unbound. In Fig. 4.3, the thresholds for binding in the dipolar
model are shown by arrows. We find that the tetramer threshold is located at
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}\approx 10$ ($h/r_{0}\approx
1.1$) and the trimer threshold, corresponding to the atom-dimer zero crossing,
occurs in the regime where all relevant length scales (scattering lengths,
dimer sizes, interaction ranges) are comparable to one another;
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}\approx 2$ ($h/r_{0}\approx
0.8$) for the dipolar model. The positions of the threshold and differences
between the results of the two models are better visible in Fig. 4.4 where we
plot the cluster energies in units of
$\hbar^{2}/ma_{{\text{B}}{\text{B}}}^{2}$.

Our numerical calculations for larger clusters indicate that, depending on
whether they are balanced ($M=N$) or not, their unbinding thresholds coincide,
respectively, with the tetramer or with the trimer ones. To understand these
results note that close to these thresholds the clusters are much larger than
the dimer. Treating there the latter as an elementary boson D, the AABBB
pentamer and the ABBB tetramer can be thought of as weakly bound DDB or DBB
“trimers” characterized by a large $a_{{\rm D}{\rm B}}$ value and repulsive DD
and BB interactions (the DD interaction is repulsive since we are above the
tetramer AABB threshold). In the limit $a_{{\rm D}{\rm B}}\rightarrow\infty$
the DD and BB interactions can be neglected and the binding energies of the
DDB and DBB composite trimers are asymptotically fractions of
$E_{{\text{A}}{\text{B}}{\text{B}}}-E_{{\text{A}}{\text{B}}}$ [99]. The ABB
trimer, ABBB tetramer, and AABBB pentamer thresholds are therefore the same
[see Fig. 4.4 (a,b,c)]. In the same reasoning, close to the AABB tetramer
crossing, the hexamer AAABBB is a weakly-bound DDD state which splits into
three dimers when the dimer-dimer attraction changes to repulsion resulting in
the same threshold value.

####  4.4.2 Threshold Determination

In this section, we numerically determine the threshold values of the few-body
clusters in the bilayer setup. To do this we need to know how the energy
depends on the interaction potential close to the threshold for unbinding. To
find out this energy dependency, let us review the principal properties of the
two-body bound state in one 1D, two 2D, and three dimensions 3D.

According to Quantum mechanics, a symmetric attractive well in 3D supports a
bound state of two particles only if the potential well depth $V$ is larger
than a critical depth $V_{c}$ [104]. Thus there is a threshold for a two-body
bound state in 3D. This is in contrast with the 1D and 2D cases where the
dimer state is formed even for infinitely small attraction between the two
particles. Therefore, in 1D and 2D the threshold for the formation of the two-
body bound state is absent [104]. In Table 4.2 we present a summary of the
principal properties of the dimer state in 1D, 2D, and 3D [105]. We notice
that for a symmetric attractive well in 2D the dimer state is weakly bound,
with its energy depending exponentially on the shallow potential $-V$,
according to

$E_{B}\approx E_{R}e^{-2cE_{R}/V},\quad\mathrm{with}\quad
E_{R}=\frac{\hbar^{2}}{mR^{2}},$ (4.11)

with c on the order of 1. This is the energy dependency we were looking for.

| 1D | 2D | 3D
---|---|---|---
$V$ | $\ll E_{R}$ | $\ll E_{R}$ | $>V_{c}\approx E_{R}$
$\psi(r>R)$ | $e^{-r/r_{B}}$ | $K_{0}(\frac{r}{r_{B}})\begin{cases}-\mathrm{log}(\frac{r}{r_{B}}),&\\!R\ll r\ll r_{B}\\\ e^{-r/r_{B}},&\\!r\gg r_{B}\end{cases}$ | $\frac{e^{-r/r_{B}}}{r}$
$r_{B}$ | $R\frac{E_{R}}{V}$ | $Re^{cE_{R}/V}$ | $R\frac{E_{R}}{V-V_{c}}$
$E_{B}=-\frac{\hbar^{2}}{mr_{B}^{2}}$ | $-\frac{V^{2}}{E_{R}}$ | $-E_{R}e^{-2cE_{R}/V}$ | $-\frac{(V-V_{c})^{2}}{E_{R}}$
Table 4.2: Bound-states in 1D, 2D and 3D for a potential well of size $R$ and
depth $V$. $\psi(r>R)$ is the wave function outside the well, $r_{B}$ is the
size of the bound state, and $E_{B}$ its energy ($E_{R}=\hbar^{2}/(mR^{2}$))
[105].

Although Eq. (4.11) is for two-body bound states we are going to use it for
larger clusters and let see if it works. Using the above result we propose to
fit the DMC binding energies with the function

$E_{{\rm A}_{N}{\rm B}_{M}}-NE_{\rm
AB}=E_{1}\exp\left[\frac{-1}{c_{1}(a_{{\text{A}}{\text{B}}}-a_{{\text{A}}{\text{B}}}^{c})+c_{2}(a_{{\text{A}}{\text{B}}}-a_{{\text{A}}{\text{B}}}^{c})^{2}}\right],$
(4.12)

for $a_{{\text{A}}{\text{B}}}>a_{{\text{A}}{\text{B}}}^{c}$, where
$a_{{\text{A}}{\text{B}}}^{c},E_{1},c_{1},c_{2}$ are free parameters. The Eq.
(4.12) can be rewritten as

$-\frac{1}{\mathrm{ln}|(E_{{\rm A}_{N}{\rm B}_{M}}-NE_{\rm
AB})/E_{1}|}=c_{1}(a_{{\text{A}}{\text{B}}}-a_{{\text{A}}{\text{B}}}^{c})+c_{2}(a_{{\text{A}}{\text{B}}}-a_{{\text{A}}{\text{B}}}^{c})^{2},$
(4.13)

which is more convenient to fit the energies. In Fig. 4.5 we show the
numerical threshold determination for the dipolar clusters. On the left panel
of Fig. 4.5 we show the threshold fitting for the trimer ABB, tetramer ABBB,
and pentamer AABBB. On the right panel of Fig. 4.5 we show the tetramer AABB
and hexamer AAABBB thresholds. The threshold values and the fitting parameters
are reported in Table 4.3. Our numerical results are consistent with our
conclusions of the previous secction, the bilayer setup have two thresholds,
one for population-imbalanced clusters at
$a^{c}_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}\approx 2$ and the
second one for population-balanced cluster at
$a^{c}_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}\approx 10$.

Figure 4.5: Fitting procedure for the determination of the threshold values of the ANBM clusters with dipolar interactions. The dashed curves correspond to the fit $c_{1}(\frac{a_{\mathrm{AB}}}{a_{\mathrm{BB}}}-\frac{a_{\mathrm{AB}}^{c}}{a_{\mathrm{BB}}})+c_{2}(\frac{a_{\mathrm{AB}}}{a_{\mathrm{BB}}}-\frac{a_{\mathrm{AB}}^{c}}{a_{\mathrm{BB}}})^{2}$. The fitting parameters are reported in Table 4.3. ANBM | $E_{1}$ | $c_{1}$ | $c_{2}$ | $a_{{\text{A}}{\text{B}}}^{c}/a_{{\text{B}}{\text{B}}}$
---|---|---|---|---
ABB | 0.0004 | 0.42307 | -0.0569969 | 2.1 $\pm$ 0.01
ABBB | 0.00041 | 0.461707 | -0.0665619 | 2.1 $\pm$ 0.05
AABBB | 0.00046 | 0.437039 | -0.042321 | 2.0 $\pm$ 0.1
AABB | 0.000089 | 0.125298 | -0.00547874 | 9.4 $\pm$ 0.2
AAABBB | 0.00026 | 0.170496 | -0.00803067 | 9.6 $\pm$ 0.2
Table 4.3: Fitting parameters for the threshold determination.

####  4.4.3 Three-Dimer Repulsion

In the last part of Subsection 4.4.1, we have integrated out the internal
degrees of freedom of the dimers, replacing them by elementary point-like
bosons. In fact, the DD zero crossing that we observe for
$a_{{\text{A}}{\text{B}}}\approx 10a_{{\text{B}}{\text{B}}}$ is a
nonperturbative phenomenon resulting from a competition between strong
repulsive and attractive interatomic forces among four individual atoms. These
interactions are strong since the corresponding scattering lengths are
comparable to the typical atomic de Broglie wave lengths $\sim 1/a_{\rm AB}$.
We emphasize that this cancellation is achieved only for two dimers. For three
dimers it is incomplete and there is a residual effective three-dimer force of
range $\sim a_{\rm AB}$ (distance, where the dimers start touching one
another). In the many-body problem, this higher-order force may compete with
the dimer-dimer interaction (if it is not completely zero) or even become
dominant. In principle, one can also discuss higher-order effects of this type
at the DB zero crossing in a DB mixture, but they are expected to be
subleading since the DD and BB interactions remain finite. In the remainder of
this section we thus concentrate on the population-balanced case.

In order to characterize the effective three-dimer interaction, we follow the
method developed previously in one dimension (Chapter 3). Namely, we analyze
the behavior of the hexamer energy just below the tetramer threshold. If the
tetramer binding energy

$E_{\rm DD}=E_{\rm AABB}-2E_{\rm AB},$ (4.14)

is much smaller than $E_{{\text{A}}{\text{B}}}$, the dimer-dimer interaction
can be considered point-like and the relative DD wave function can be
approximated by

$\phi(r)\propto K_{0}(\kappa r),$ (4.15)

where $\kappa=\sqrt{-2mE_{\rm DD}/\hbar^{2}}$ is the inverse size of the
tetramer. Similarly, the AAABBB hexamer under these conditions reduces to the
well-studied problem of three point-like bosons [106, 107, 108, 109, 110, 111,
98], according to which the ground-state hexamer binding energy

$E_{\rm DDD}=E_{\rm AAABBB}-3E_{\rm AB},$ (4.16)

should satisfy [110, 111]

$E_{\rm DDD}/|E_{\rm DD}|=-16.5226874.$ (4.17)

We expect the ratio $E_{\rm DDD}/|E_{\rm DD}|$ to reach the zero-range limit
(4.17) as we approach the dimer-dimer zero crossing, i.e., as $\kappa a_{\rm
AB}\rightarrow 0$. In Fig. 4.6, we plot $E_{\rm DDD}/|E_{\rm DD}|$ versus
$\kappa a_{\rm AB}$ and indeed see a tendency towards the value (4.17)
although the effects of the finite size of the dimers and their internal
degrees of freedom, that we have neglected in the zero-range model, are
obviously important. The fact that the hexamer energy lies above the limit
(4.17) points to an effective three-dimer repulsive force. We note again that
the values of the ratio $E_{\rm DDD}/|E_{\rm DD}|$ obtained for Gaussian and
dipolar potentials are quite close to each other for all values of
$a_{{\text{A}}{\text{B}}}$ suggesting a certain universality of this problem
and a relative unimportance of the long-range interaction tails.

In order to quantify the three-dimer interaction observed in Fig. 4.6, we
compared our results with a zero-range model with hard-core hyperradius
constraint $\rho_{0}$, developed in Ref. [71].

Figure 4.6: The hexamer-to-tetramer binding energy ratio $E_{\rm DDD}/|E_{\rm
DD}|$ as a function of $\kappa a_{\rm AB}$. The solid line is the result of
the zero-range model with the hard-core hyperradius constraint at
$\rho_{0}=0.13a_{\rm AB}$ [71].

The authors in Ref. [71] extended the model of three point-like dimers by
requiring that the three-dimer wave function vanishes at a hyperradius
$\rho_{0}$. For three dimers, with coordinates ${\bf r}_{1}$, ${\bf r}_{2}$,
and ${\bf r}_{3}$, the hyperradius is defined as

$\rho=\sqrt{x^{2}+y^{2}},$ (4.18)

where

${\bf x}=\frac{2{\bf r}_{3}-{\bf r}_{1}-{\bf
r}_{2}}{\sqrt{3}}\quad\mathrm{and}\quad{\bf y}={\bf r}_{1}-{\bf r}_{2},$
(4.19)

are the Jacobi coordinates. For this minimalistic model $E_{\rm DDD}/|E_{\rm
DD}|$ is a function of the ratio $\kappa\rho_{0}$, relating the three- and
two-dimer interaction strengths. Kartavtsev and Malykh [111] discussed the
adiabatic hyperspherical method in detail and applied it to the $\rho_{0}=0$
limit, i.e., the case of no three-body interaction. The authors in Ref. [71]
accounted for finite $\rho_{0}$, by set the hyperradial channel functions to
zero at $\rho=\rho_{0}$. In this way, they obtained the ratio $E_{\rm
DDD}/|E_{\rm DD}|$ as a function of $\kappa\rho_{0}$ (for more details see
Ref. [71]). We then treat $\rho_{0}$ as a constant (independent of $E_{\rm
DD}$) determined by fitting the DMC and SVM data in $\kappa a_{\rm AB}<0.4$
range. By minimizing $\chi^{2}$ we obtain $\rho_{0}=0.13a_{\rm AB}$.

The inclusion of the three-body hard-core constraint, even corresponding to
numerically very small $\kappa\rho_{0}$, leads to a spectacular deviation from
Eq. (4.17). This interesting effect is due to an enhancement of the three-
dimer interaction by strong two-dimer correlations.

Promising candidates for observing the predicted cluster states are bosonic
dipolar molecules characterized by large and tunable dipolar lengths which, at
large electric fields, tend to $r_{0}=5\times 10^{-6}$m for 87Rb133Cs [112,
113], $r_{0}=2\times 10^{-5}$m for 23Na87Rb [114, 115] and $r_{0}=6\times
10^{-5}$m for 7Li133Cs [116]. Fermionic 87Rb40K [117, 118] and 23Na40K [119,
120, 121] molecules ($r_{0}=7\times 10^{-7}$m and $r_{0}=7\times 10^{-6}$,
respectively) could be turned into bosons by choosing another isotope of K.
The interlayer distance, fixed by the laser wavelength, has typical values of
$h\approx(2-5)\times 10^{-7}$m, which is thus sufficient for observing the
few-body bound states that we predict for ratios $h/r_{0}>0.8$.

The next step in this work is to generalize these findings to the many-body
problem when a new scale (density $n$) comes into play. It is important to
understand how the two- and three-body effects correlate with each other as
one passes through the dimer-dimer zero crossing. Although we find no
qualitative difference between the dipolar and Gaussian models in our few-body
results, the long-range tails will be important when the quantity $nr_{0}$
becomes comparable to the inverse healing length (which is where the dipolar
condensate becomes rotonized). For bilayer dipoles the relevant region of
parameters is close to the dimer-dimer zero crossing, which we predict to be
at $h/r_{0}\approx 1.1$. In Chapter 6, we studied the many-body problem of
dipolar bosons in a bilayer geometry.

###  4.5 Summary

To summarize, we have studied few-body clusters ANBM with $N\leq M\leq 3$ in a
two-dimensional Bose-Bose mixture using different (long-range dipolar and
short-range Gaussian) intraspecies repulsion and interspecies attraction
models. In both cases, the intraspecies scattering length
$a_{{\text{A}}{\text{A}}}=a_{{\text{B}}{\text{B}}}$ is of the order of the
potential ranges, whereas we tune $a_{{\text{A}}{\text{B}}}$ by adjusting the
AB attractive potential (or the interlayer distance in the bilayer setup). We
find that for $a_{{\text{A}}{\text{B}}}\gg a_{{\text{B}}{\text{B}}}$ all
considered clusters are (weakly) bound and their energies are independent of
the interaction model. As the ratio
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}$ decreases, the increasing
intraspecies repulsion pushes the clusters upwards in energy and eventually
breaks them up into $N$ dimers and $M-N$ free B atoms. In the population
balanced case ($N=M$) this happens at
$a_{{\text{A}}{\text{B}}}/a_{{\text{B}}{\text{B}}}\approx 10$ where the dimer-
dimer attraction changes to repulsion. By studying the AAABBB hexamer near the
dimer-dimer zero crossing we find that it very much behaves like a system of
three point-like particles (dimers) characterized by an effective three-dimer
repulsion. A dipolar system in a bilayer geometry can thus exhibit the
tunability and mechanical stability necessary for observing dilute liquids and
supersolid phases.

##  Chapter 5 Quantum Halo States in Two-Dimensional Dipolar Clusters

The purpose of the present chapter is to study the ground-state properties of
loosely bound dipolar clusters composed of two to six particles in a two-
dimensional bilayer geometry. We investigate whether halos, bound states with
a wave function that extends deeply into the classically forbidden region
[122, 123], can occur in this system. The dipoles are confined to two layers,
A and B, with dipolar moments aligned perpendicularly to the planes. The
binding energies, pair correlation functions, spatial distributions, and sizes
are calculated at different values of the interlayer distance by using the
diffusion Monte Carlo method. We find that for large interlayer separations
the AB dimers are halo states and following a universal scaling law relating
the energy and size of the bound state. For ABB trimers and AABB tetramers, we
find two very distinct halo structures. For large values of the interlayer
separation, such halo states are weakly bound and the typical distances
between BB and AB dipoles are similar. However, for the deepest bound ABB and
AABB clusters, and as the clusters approach the unbinding energy threshold, we
find a highly anisotropic structure, in which the AB pair is on average
several times closer than the BB pair. Similar symmetric and asymmetric
structures are observed in pentamers and hexamers, both being halo states,
thus providing halos with the largest number of particles ever observed or
predicted before.

###  5.1 Introduction

One of the most remarkable aspects of ultracold quantum gases is their
versatility, which permits to bring ideas from other areas of physics and
implement them in a clean and highly controllable manner. Some of the examples
of fruitful interdisciplinary borrowings include Efimov, states originally
introduced in nuclear physics and observed in alkali atoms [124, 125, 126],
lattices created with counter-propagating laser beams [127, 128] as models for
crystals in condensed matter physics, Bardeen-Cooper Schrieffer (BCS) pairing
theory first introduced to explain superconductivity and later used to