##  Chapter 0 Path Integral Monte Carlo for bosons in continuous space

Monte Carlo methods have become a popular approach to tackling physics
problems. Although the statistics community was already promoting the
technique in the early 1900s, it was Metropolis1949 who, 50 years later111Of
course, it became more practical with the rise of computers., first published
the approach in a physical context to compute the ground-state energy of the
Schrödinger equation Reynolds1990222Actually, Stanislaw Ulam and John von
Neumann used the Markov Chain Monte Carlo method when working on neutron
diffusion for nuclear weapons in the Manhattan Project at the Los Alamos
National Laboratory Metropolis1987; Benov2018.. The idea was inspired by
Fermi’s remarks given at a symposium Fermi1957 where he noticed the
resemblance of the time-dependent Schrödinger equation to a diffusion equation
(see later) and hence can be treated as a stochastic process. Diffusion Monte
Carlo was discovered.

Nowadays, Monte Carlo methods are used in all fields of physics, from
molecular modelling in biophysics Carlon2002 to designing detectors in
particle physics Agostini2018 or understanding galaxy evolution in
astrophysics MacGillivray1982. In statistical physics, the overall purpose of
using the Monte Carlo method is to compute a multivariate integral over phase
space to compute the expectation value of some macroscopic variable. The
motivation being that the error converges as $\sim 1/\sqrt{n}$, where $n$
represents the total number of independent samples333In practice, the samples
are correlated (see section 5)., regardless of the dimension of the phase
space. This property is especially advantageous for quantum systems where the
Hilbert space expands exponentially with the number of particles. Other
methods, such as exact diagonalization, thereby scale very unfavourably with
system size due to memory requirements. By contrast, _Quantum Monte Carlo_
(QMC) algorithms are computational methods developed to produce solutions to
quantum many-body systems by handling the multivariate integrals that emerge
in various formulations of the many-body problem.

The key disadvantage of the QMC approach is that it is almost always only
effective for bosonic systems444There are some exceptions to this rule, e.g.,
bold diagrammatic Monte Carlo Prokofev2007.. The notorious numerical _sign
problem_ , in which near-cancellation of the positive and negative
contributions to the integral, prevents general fermionic systems from being
sampled in this way Pollet2012. This does not prohibit the Monte Carlo
approach from being used, but it does restore the exponential scaling of
computational resources with system size. Solving the fermion-sign problem is
seen to be impossible, grounded on the fact that it has been proven to be non-
deterministic polynomial (NP)-hard Troyer2005.

The QMC algorithm we will be using is called Path Integral (Quantum) Monte
Carlo, or PIMC. The method was developed by Pollock and Ceperley Ceperley1984;
Ceperley1987; Ceperley1995 to understand the $4$He superfluid transition.
Nowadays, it is used to understand a wide variety of condensed matter systems,
e.g., dipolar systems Cinti2017, ultracold gases Cinti2014; Pascual2021;
Carleo2013 and quantum fluids and solids Boninsegni2009; Corboz2008. By
developing the method, Pollock and Ceperley were able to turn Feynman’s theory
of superfluidity Feynman1953 into a quantitative theory that could be
validated using experimental evidence. However, the method was still limited
by the number of particles in the system, as the sampling of Bose statistics
becomes inefficient at low temperatures. Nevertheless, the hurdle was overcome
by using an alternative sampling approach, known as the worm algorithm (WA)
Prokofev1998; Boninsegn2006; Boninsegn2006PRL. The combination of both methods
is often called Worm Algorithm Path integral Monte Carlo (WAPIMC). In the case
of ultracold atoms in optical lattices, WAPIMC has been shown to be in good
agreement with experiment Trotzky2010. WAPIMC allows for a lot of
thermodynamic observables of interest at temperatures greater than zero,
including structural, e.g., pair correlation functions, the structure factor,
and compressible, as well as quantities such as superfluid density and
condensate fraction.

In this chapter, we will explain the WAPIMC method in continuous space,
starting with the basics of the PIMC algorithm (1) and thereby also Markov
chain Monte Carlo (3). Subsequently, in section 2, a discussion of the small-
time propagator can be found. After, a depiction of the sampling updates (3)
is given, with a special focus on the Worm Algorithm (2). Finally, a brief
overview of the computation of physical observables (4) is given, complemented
by an explanation of Monte Carlo data error analysis (5), including
autocorrelation effects.

###  1 Path Integral Monte Carlo

In the following subsections, WAPIMC is presented and follows the approach
used by Yao2020; Ceperley1995; Boninsegn2006. Some other useful didactic
sources considered are Pollet2012; YanBlume2017; dornheim2014.

####  1 Feynman’s path integral

Consider a system of $N$ identical _bosons_ in a volume $V$ of $d$ dimensions
at a finite temperature $T$, which is represented in the canonical ensemble by
a Hamiltonian $\ham$. Quantum statistical mechanics feynman1965; feynman1972
dictates that the expectation value of an observable $\hat{O}$ at
thermodynamic equilibrium is given by

$\displaystyle\expval*{\hat{O}}=\frac{1}{Z}\Tr[\hat{\rho}\
\oper]\qq{with}Z=\Tr[\hat{\rho}]\qq{and}\hat{\rho}=e^{\beta\ham}$

where $Z$ is the canonical partition function, $\hat{\rho}$ the matrix density
operator, and $\beta=(k_{B}T)^{-1}$ the inverse temperature, with $k_{B}$ the
Boltzmann constant. In order to evaluate the trace over the configurations,
one has to choose a convenient basis. Let us take the spatial locations of the
N particles
$\ket{\mathbf{R}}\equiv\ket{\mathbf{r}_{1},\cdots,\mathbf{r}_{N}}$, so that
the trace for an operator $\hat{X}$ can be written as

$\displaystyle\operatorname{Tr}[\hat{X}]=\int\dd{\mathbf{R}}\expval{\hat{X}\hat{S}}{\mathbf{R}}\qq{where}\hat{S}=\frac{1}{N!}\sum_{\sigma\in\mathcal{P}}\ket{\sigma\cdot\mathbf{R}}\bra{\mathbf{R}}$

The symmetrization operator $\hat{S}$ must be included in the trace to account
for the indistinguishability of the identical bosons and, thereby, a crucial
part of the physics of the system. The canonical partition function $Z$ is
then written by

$\displaystyle
Z=\frac{1}{N!}\sum_{\sigma\in\Pi}\int\dd{\mathbf{R}}\mel{\mathbf{R}}{\
\hat{\rho}\
}{\sigma\cdot\mathbf{R}}=\frac{1}{N!}\sum_{\sigma\in\mathcal{P}}\int\dd{\mathbf{R}}\rho(\mathbf{R},\sigma\cdot\mathbf{R};\beta)$

where $\rho(\mathbf{R},\mathbf{R}^{\prime};\beta)$ is called the _density
matrix_ at inverse temperature $\beta$. If the Hamiltonian $\ham$ is
hermitian, the density matrix is symmetric in its first two arguments.

Often, the potential energy555My apologies for the unconventional notation.
Nevertheless, this part of the contribution to the density matrix is often
referred to as the action (see section 2). In the Hamiltonian, the potential
for the optical lattice and the interaction will be noted as V and U,
respectively. part $\hat{\mathcal{S}}$ of the Hamiltonian
$\ham=\hat{\mathcal{T}}+\hat{\mathcal{S}}$ is diagonal with respect to the
position eigenbasis. Hence, it would be preferable to decompose the density
matrix into a product containing the kinetic $\hat{\mathcal{T}}$ and potential
$\hat{\mathcal{S}}$ contributions to the Hamiltonian. Nevertheless, in
general, we have that these contributions do not commute. One could remedy
this problem by using the Baker-Campbell-Hausdorff (BCH) formula Campbell1897
to construct an approximation given by

$\displaystyle
e^{-\beta\ham}=e^{-\beta\hat{\mathcal{T}}}e^{-\beta\hat{\mathcal{S}}}e^{-\frac{\beta^{2}}{2}\comm{\hat{\mathcal{T}}}{\hat{\mathcal{S}}}\
+\
\cdots}=e^{-\beta\hat{\mathcal{T}}}e^{-\beta\hat{\mathcal{S}}}+\mathcal{O}(\beta^{2})$

Splitting the kinetic and potential contribution with the help of BCH formula
is often called the _primitive approximation_. Notice that, the error to our
approximation diverges in the low temperature (quantum) limit $\beta\gg 1$,
the limit one is often interested in. To remedy this problem, one can split
the density matrix into $J\in\mathbb{Z}$ parts of step size
$\tau\equiv\beta/J$ by using the normalization condition
$\hat{I}=\int\dd{\mathbf{R}}\ket{\mathbf{R}}\bra{\mathbf{R}}$. Doing this
yields

$\displaystyle
Z=\frac{1}{N!}\sum_{\sigma\in\Pi}\int\dd{\mathbf{R}_{J-1}}\mydots\dd{\mathbf{R}_{0}}\mel{\mathbf{R}_{0}}{e^{-\tau\hat{H}}}{\mathbf{R}_{1}}\mel{\mathbf{R}_{1}}{e^{-\tau\hat{H}}}{\mathbf{R}_{2}}\mydots\mel{\mathbf{R}_{J-1}}{e^{-\tau\hat{H}}}{\sigma\cdot\mathbf{R}_{0}},$
(1)

where we introduce the notation666In general, I would like to reserve the
indices $i$ and $j$ for enumerating the number of particles and time slices
respectively. The index $k$ can be used as a spare index.
$\ket{\mathbf{R}_{j}}\equiv\ket{\mathbf{r}_{1,j},\cdots,\mathbf{r}_{N,j}}$ and
$\sigma$ is a permutation operator of the group of permutations of the $N$
elements $\Pi$. Notice that for the above expression, we used that the
symmetrization operator $\hat{S}$ commutes with the Hamiltonian $\ham$.
Furthermore, we define
$\ket{\mathbf{R}_{J}}\equiv\ket{\sigma\cdot\mathbf{R}_{0}}$ such that we can
write the partition function in a more compact form:

$\displaystyle
Z=\frac{1}{N!}\sum_{\sigma\in\Pi}\prod_{j=0}^{J-1}\int\dd{\mathbf{R}_{j}}\rho(\mathbf{R}_{j},\mathbf{R}_{j+1};\tau).$
(2)

The expression for the partition function is still exact for every $J\geq 1$.
Although, we modified the problem from integrating one density matrix at
temperature $\beta$ to integrating $J$ density matrices at $\tau$. Now it will
be easier to find an expression for the density matrix. In practice, one takes
$J$ very large such that $\tau\ll 1$ and an accurate expression of
$\rho(\mathbf{R}_{j},\mathbf{R}_{j+1};\tau)$ can be found using (1) and
higher-orders (see section 2). A broader discussion is given in section 2.

A key insight into the problem can be made if one identifies $\beta=it/\hbar$
such that the operator $\operatorname{exp}(-\beta\ham)$ turns into the unitary
time evolution operator. Consequently, the system can be seen as N particles
evolving through imaginary time, i.e. the density matrix becomes a _time
evolution propagator_ 777The same substitution can be made in the time-
dependent Schrödinger equation which results in a diffusion equation with a
diffusion constant $D=\tfrac{\hbar^{2}}{2m}$.. This mapping of a
$d$-dimensional quantum system to a $d+1$ classical system is known as the
quantum-classical mapping of Feynman path integral formulation of quantum
statistical mechanics. From now on, we shall use both representations
interchangeably.

Figure 1: Two representations of the worldline configurations. The schematic
are inspired by and adapted from Yao2020; dornheim2014. red(a) A collection of
beads connected by links plotted in function of imaginary time. red(b) Two-
dimensional spatial position of the beads connected by their links to form a
ring or polymer.

A configuration $\config$ is defined as the collection
$(\mathbf{R}_{1},\cdots,\mathbf{R}_{J})$ which represents $N$ discrete
worldlines in imaginary time from 0 to $\beta$. An example of a typical
configuration in one spatial dimension can be found in figure 1red(a). The
position of the i-th particle at the j-th time slice, i.e. $r_{i,j}$ is called
a _bead_ and the propagator
$\rho\left(\mathbf{r}_{i,j},\mathbf{r}_{i,j+1};\tau\right)$ who connects two
beads of consecutive time slices is called a _link_. In table 1, a list of all
the jargon of a PIMC simulation van be found. As can be seen in (1), due to
the imaginary time boundary condition
$\ket{\mathbf{R}_{J}}=\ket{\sigma\cdot\mathbf{R}_{0}}$, the worldlines of are
$n\beta$ periodic with $1\leq n\leq N$. In the case $n>1$, the particles are
interchanged in the configuration, which accounts for the indistinguishability
of the particles.

ChandlerWolynes1981 introduced another perspective. In the classical
representation, one could view the worldlines as ring polymers where they form
a polyatomic fluid, see figure 1red(b). As the particle permutes, the polymers
grow with a maximum of $NJ$ beads. It has been shown that the polymer spatial
area is proportional to the de Broglie wavelength Ceperley1995

$\displaystyle\lambda_{\beta}=\frac{2\pi\hbar}{\sqrt{2\pi
mk_{\mathrm{B}}T}}=\sqrt{4\pi\lambda\beta},$ (3)

with $\lambda=\hbar^{2}/2m$. Hence, the quantum statistics become important
when the thermal de Broglie wavelength of the particles is of the same order
as the inter-particle distance of the particles and, thereby, the polymers can
permute, i.e. the particles become indistinguishable. Accordingly, if the
temperature is raised, i.e., $\beta\to 0$, the polymer shrinks and becomes
distinguishable again, analogously to the single-particle quantum
wavefunction.

With the procedure of constructing (1) in mind, equation (1) can be written in
an insightful form. Namely, one can define the normalized _probability density
function_ (pdf) of a configuration $\config$ by

$\displaystyle\pi(\config)\equiv\frac{\mathcal{W}(\config)}{Z}=\frac{1}{Z}\frac{1}{N!}\mel{\mathbf{R}_{J}}{e^{-\tau\hat{H}}}{\mathbf{R}_{J-1}}\ldots\mel{\mathbf{R}_{1}}{e^{-\tau\hat{H}}}{\mathbf{R}_{0}},$
(4)

where it satisfies the normalization condition
$\int\dd{\config}\pi(\config)=1$ and $\mathcal{W}(\config)$ is called the
weight of a configuration $\config$. Hence, we write the expectation value of
the observable $\oper$ as

$\displaystyle\expval*{\hat{O}}=\int\dd{\config}\pi(\config)\mathcal{Q}(\config)\qq{with}\mathcal{Q}(\config)=\frac{\mel{\mathbf{R}_{1}}{e^{-\tau\hat{H}}\oper}{\mathbf{R}_{0}}}{\mel{\mathbf{R}_{1}}{e^{-\tau\hat{H}}}{\mathbf{R}_{0}}},$
(5)

where $\mathcal{Q}(\config)$ is called the _estimator_ for observable $\oper$
associated with the configuration $\config$. The pdf $\pi(\config)$ and
estimator $\mathcal{Q}(\config)$ are worked out in sections 2 and 4,
respectively. Nevertheless, evaluating these identities for all configurations
would be unfeasible. This can be solved using Monte Carlo integration.

####  2 Monte Carlo integration

Consider a function $f$ from $\mathbb{R}^{d}$ to $\mathbb{R}$ where
$d\in\mathbb{Z}$ and the integral

$\displaystyle I=\int_{\mathbb{R}^{d}}f(\vb{x})\dd{\vb{x}}$

The integral is often difficult to compute in high-dimensional integration
spaces. A common way is to divide the integration space into small cells, each
with its own evaluation point $\vb{x}_{i}$, and then sum the value of the
function at these locations multiplied by thee volume of the cells
$\dd{\vb{x}_{i}}$, i.e. $I\approx\sum_{i}f(\vb{x}_{i})\dd{\vb{x}_{i}}$.
Nevertheless, the computational cost scales exponentially with the dimension
$d$. For equation 5 the dimension of integration space, i.e., the
configuration space, $D=dNJ$ becomes very large if one wants to find an
accurate small-time propagator $\rho(\mathbf{R},\mathbf{R}^{\prime};\tau)$ and
thus taking $J\gg 1$.

One of the methods to prevail the so-called curse of dimensionality is Monte
Carlo integration, which samples the configuration stochastically. The key
reason for the above method’s failure is that it treats every part of the
configuration on the same footing. Conversely, Monte Carlo integration uses
importance sampling, where the configuration with a large contribution has a
higher chance of being generated than the samples with a low contribution.
More concretely, suppose the stochastic variable $X$ which lies in
$\mathbb{R^{D}}$ distributed by a pdf $\pi$, one writes

$\displaystyle
I=\int_{\mathbb{R}^{D}}\frac{f(x)}{\pi(x)}\pi(x)\dd{x}=\langle\mathcal{I}(X)\rangle_{\pi}\approx\frac{\mathcal{I}\left(X_{1}\right)+\mathcal{I}\left(X_{2}\right)+\ldots+\mathcal{I}\left(X_{n}\right)}{n},$
(6)

where $\mathcal{I}(X)=f(X)/\pi(X)$. The outcome of the integral is thus the
expectation value of $\mathcal{I}(x)$ under the pdf $\pi$. Consequently,
because of the law of large numbers, the integral can be approximated by a
large set of independent samples $X_{1},X_{2},\ldots,X_{n}$ as shown in (6).
The central limit theorem tells us that the error will be,
$\varepsilon_{I}=\sigma_{I}/\sqrt{n}$ where $\sigma_{I}$ is the standard
deviation of $\mathcal{I}(x)$. In section 5 a more in-depth analysis of the
error of PIMC will be given.

If one makes the identification with equation 5 obtained in the previous
section, we find that the expectation value can be approximated by

$\displaystyle\expval*{\oper}\approx\frac{\mathcal{Q}\left(\config_{1}\right)+\ldots+\mathcal{Q}\left(\config_{n}\right)}{n},$
(7)

where $\config_{1},\ldots,\config_{n}$ are sampled out of pdf $\pi$ given by
(4). However, the direct sampling of the configuration $\config$ out of a pdf
as complex as $\pi$ can be a hard task. The solution is sampled with the help
of a Markov chain, which will be discussed in the next section.

####  3 Markov Chain Monte Carlo

A _Markov chain_ is a discrete stochastic sequence, in our case, of
configurations that starts from an initial configuration $\config$ and
generates new configurations with the help of a set of updates, each with
their own transition probability $\mathcal{T}(\config\to\config^{\prime})$.
The result is a chain of configurations

$\displaystyle\config_{1}\rightarrow\config_{2}\rightarrow\ldots\rightarrow\config_{k}\rightarrow\ldots\
.$ (8)

Each transition is independent of the history of the chain, i.e., they satisfy
_Markov’s property_ , and the normalization condition
$\int\dd{\config^{\prime}}\mathcal{T}(\config\to\config^{\prime})=1$. The goal
is to generate a Markov chain where the _equilibrium distribution_ of the
configurations in the chain is the pdf $\pi$ given by equation 4. To achieve
this, the Markov chain performs a random walk in the configuration space that
is attracted to configurations with high contributions to the integral (5).

For the walker to produce a stationary distribution, one usually requires two
conditions: _detailed balance_ and _ergodicity_. The former states that the
probability of being in a configuration $\config$ times the transition
probability of $\config\to\config^{\prime}$ is equal to the probability of
being in a configuration $\config^{\prime}$ times the transition probability
of $\config^{\prime}\to\config$, i.e.,

$\displaystyle\pi(\config)\mathcal{T}\left(\config\rightarrow\config^{\prime}\right)=\pi\left(\config^{\prime}\right)\mathcal{T}\left(\config^{\prime}\rightarrow\config\right).$
(9)

In other words, each transition is reversible. The latter condition states
that given two arbitrary configurations $\config$ and $\config^{\prime}$ there
always exists a path or sequence of transitions with a combined non-zero
transition probability connecting them. It ensures that the whole
configuration space can be sampled. Note that detailed balance and ergodicity
together are sufficient for the Markov chain to generate a stationary
distribution, but not necessary. There are also walkers which do not satisfy
detail balance and produce stationary distributions that correspond with the
desired distribution $\pi$. A practical example is given in Bernard2009.

The walker of interest for our purposes is the _Metropolis-Hastings
algorithm_. It is constructed in such a way that it satisfies the above two
conditions, and follows two steps:

  1. [(i)]

  2. 1.

First, propose a move or update which modifies the last element of the Markov
chain such that a new configuration $\config^{\prime}$ is offered with
_proposal probability_
$\mathcal{P}\left(\config_{i}\rightarrow\config^{\prime}\right)$. How this
update is constructed is discussed in section 3.

  3. 2.

Afterwards, the new configuration $\config^{\prime}$ is evaluated to check
whether we want to accept or reject the move, determined by the _acceptance
probability_

$\mathcal{A}\left(\config_{i}\rightarrow\config^{\prime}\right)=\min\left(1,\frac{\pi\left(\config^{\prime}\right)\mathcal{P}\left(\config^{\prime}\rightarrow\config_{i}\right)}{\pi(\config_{i})\mathcal{P}\left(\config_{i}\rightarrow\config^{\prime}\right)}\right),$
(10)

such that the _transition probability_
$\mathcal{T}(\config_{t}\to\config^{\prime})=\mathcal{P}\left(\config_{i}\rightarrow\config^{\prime}\right)\mathcal{A}\left(\config_{i}\rightarrow\config^{\prime}\right)$.
If it is accepted, one sets $\config_{i+1}=\config^{\prime}$, otherwise,
$\config_{i+1}=\config_{i}$.

Together, the steps 1 and 2 lead to a new configuration. The acceptance
probability (10) can easily be checked to satisfy the detailed balance
equation by substituting it into equation 9 and considering the two cases
$\pi\left(\config^{\prime}\right)\mathcal{P}\left(\config^{\prime}\rightarrow\config_{i}\right)<\pi(\config_{i})\mathcal{P}\left(\config_{i}\rightarrow\config^{\prime}\right)$
and
$\pi(\config_{i})\mathcal{P}\left(\config_{i}\rightarrow\config^{\prime}\right)<\pi\left(\config^{\prime}\right)\mathcal{P}\left(\config^{\prime}\rightarrow\config_{i}\right)$888Note
that the partition function $Z$ is the same for both configuration
$\config^{\prime}$ and $\config_{i}$, and hence, thereby cancels..

####  4 PIMC Algorithm

In summary, a common PIMC simulation is formulated as follows.

  1. [(i)]

  2. 1.

By using the detailed balance equation (9) an ergodic set of updates can be
derived, each with its own proposal probability $\mathcal{P}$ (see section 3).
Note that the ergodicity of the set is system dependent.

  3. 2.

Compute an (approximated) expression for the density matrix with the help of
the Hamiltonian $\ham$ which yields an exact result in the limit $\tau\to 0$
(see section 2).

  4. 3.

Start from an arbitrary initial configuration $\config_{k}$ with $k=0$ that
respect the system parameters set by the Hamiltonian $\ham$.

  5. 4.

Create a new configuration $\config^{\prime}$ with a random update from the
set constructed in (i) and evaluate the acceptance probability
$\mathcal{A}\left(\config_{k}\rightarrow\config^{\prime}\right)$ given by
equation 10.

  6. 5.

If the acceptance probability
$\mathcal{A}\left(\config_{k}\rightarrow\config^{\prime}\right)$ is larger
than a random number $a\in[0,1)$, the new configuration is added to the chain,
i.e., $\config_{k+1}=\config^{\prime}$, else the initial configuration is
added to the chain, or, $\config_{k+1}=\config_{k}$.

  7. 6.

Last, the estimators $\mathcal{Q}$ of the desired collection of observables
$\\{\hat{O}_{i}\\}$ (see section 4), defined in (5), are evaluated for the
added member of the chain $\config_{k+1}$.

  8. 7.

Set $k$ to $k+1$.

This process of steps (iii-vii) is then repeated until enough Markov chain
entries are collected. In the case of a rejection, the old configuration must
be still counted as a member of the chain. Including only accepted
configurations leads to wrong results.

Table 1: Jargon regarding the Path integral Monte Carlo simulation. [()]15pt Bead | $\mathbf{r}_{k,j}$ | a single coordinate of the $k$-th particle at the $j$-th imaginary time index
---|---|---
[()]20pt Time slice | $\mathbf{R}_{j}$ | a set of beads at the $j$-th imaginary time index
[()]20pt Configuration | $\left\\{\mathbf{R}_{0},\cdots,\mathbf{R}_{J}\right\\}$ | the set of all time slices which as a whole defines the configuration
[()]20pt Link | $\rho\left(\mathbf{R}_{j},\mathbf{R}_{j+1};\tau\right)$ | the propagator connecting two consecutive time slices
[()]20pt Single-particle link | $\rho\left(\mathbf{r}_{k,j},\mathbf{r}_{k,j+1};\tau\right)$ | the propagator connecting two consecutive beads
[()]20pt Weight function | $\mathcal{W}(\config)$ | the weight of configuration $\config$
[()]20pt Probability density function | $\pi(\config)$ | the normalized probability to be in configurations $\config$
[()]20pt Estimator function | $\mathcal{Q}(\config)$ | the estimator of observable to be evaluated for each configuration $\config$ in the Markov chain
[()]20pt Transition probability | $\mathcal{T}\left(\config\rightarrow\config^{\prime}\right)$ | the probability a configuration $\config$ transitions to configuration $\config^{\prime}$
[()]20pt Proposal probability | $\mathcal{P}\left(\config\rightarrow\config^{\prime}\right)$ | the probability a configuration $\config^{\prime}$ is proposed when the system is in $\config$
[()]20pt Acceptance probability | $\mathcal{A}\left(\config\rightarrow\config^{\prime}\right)$ | the probability a configuration $\config^{\prime}$ is accepted when the system is in $\config$

##### Comparison with other methods

The last part of this section I would like to devote to how PIMC999We use PIMC
to refer to the core method of the simulation, as WA only adds additional
updates. compares to other methods and why it is the method of choice for our
purposes. In doing this, Ceperley1995; YanBlume2017; king_scaling_2021;
Pollet2012 were considered.

First and foremost, PIMC is used by a lot of others to simulate these
ultracold atoms in optical lattices Cinti2014; Pascual2021; Carleo2013;
Gauter2021101010Which is a priori, of course, not a strong argument on why to
use PIMC.. It is well accepted as the standard in the field for computing the
accurate estimate of finite-temperature equilibrium statistics with dimensions
higher than one king_scaling_2021. That said, the main reason why it has this
status is that it is completely unbiased. The only input PIMC needs is a
microscopic Hamiltonian and does not require any a priori assumptions, i.e.,
it is an ab-initio simulation. By contrast, comparable methods such as Path
Integral Ground State (PIGS) YanBlume2017, Variational Monte Carlo (VMC),
Green’s-function Monte Carlo (GFMC) and Diffusion Monte Carlo (DMC) need a
trail wavefunction to start from. However, in general, when humans appreciate
the outcome of a simulation, they tend to not look deeper into the result.
Hence, a systematic inaccuracy is unavoidably introduced in the above-
mentioned methods. Moreover, for complicated systems, a good quality trail
function has far too many potential parameters, making it difficult to
optimize properly. One may refer to PIMC as _numerical exact_ for boson
systems, which means that all the approximations required for the fact that it
is a numerical technique can be decreased arbitrarily as long as sufficient
computing is performed. Although one has to remember that for fermions, one
encounters the Fermi sign problem. Another advantage of PIMC is that
observables related to particle statistics are easier expressed with
comparison to other methods Ceperley1995.

The disadvantage of PIMC is that it is quite involved compared to other
methods such as VMC, both in structure and computationally. VMC is simple to
comprehend and to program. As a result, for a system with for example a lot of
symmetry, it can be preferred. In addition, the particle statistics remain
difficult for WAPIMC. Sampling worldlines with many-particle permutations
still scales unfavourably with the number of particles (see section 2).
Although WA greatly improved the problem, it remains the main obstacle even
when one uses WA Pollet2012. Also, note that WAPIMC does not have the ability
to measure strictly zero-temperature systems. Nevertheless, in the context of
comparing with experimental data, this is no problem.

A similar formulation of PIMC exists, where one restricts the spatial location
of the particles to the vertices of the lattice. In doing this, one can
formulate the algorithm such that there is no time step error. However, in
practice, the time step error does not seem to be the main bottleneck of the
former, i.e., the more ab-initio continuous space PIMC variant.

###  2 Small-time propagator

One obstacle with working in Feynman’s path integral formulation is that one
has an expression in terms of the many-body density matrix. However, this
object is unknown in exact form for any non-trivial system. Hence, one has to
find good approximations. Depending on the system at hand, one has to
construct clever approximations which are valid in the correct regime. For
ultracold atoms in optical lattices, the Hamiltonian is written as follows:

$\displaystyle\ham=\sum_{i=1}^{N}-\lambda\hat{\nabla_{\mathrm{i}}}^{2}+\sum_{i=1}^{N}V\left(\hat{\mathbf{r}}_{i}\right)+\sum^{N-1}_{i=1}\sum^{N}_{k=i+1}U\left(\hat{\mathbf{r}}_{i}-\hat{\mathbf{r}}_{k}\right),$
(23)

where we defined $\lambda=\frac{\hbar^{2}}{2m}$ for convenience. The first
term on the right side of the equality represents the kinetic energy
contribution, V the optical landscape or external potential at a given
position, and U the interaction between the relative positions of two
particles $i$ and $k$, i.e.,
$\hat{\vb{r}}_{ik}=\hat{\mathbf{r}}_{i}-\hat{\mathbf{r}}_{k}$. As discussed in
chapter LABEL:ch:_Ultracold_Atoms_In_Optical_Lattices, ultracold gases are
often diluted gases. As the particles rarely collide, one can consider only
the two-body interaction which is done in the last term of the Hamiltonian.
This is the first approximation we make, nevertheless, a good experimentally
confirmed one Bloch2008.

Considering that the interactions only depend on the relative position of the
particles, it is favourable to extract the relative Hamiltonian of the reduced
single particle from equation (23), i.e.,
$\ham^{\mathrm{rel}}=\ham_{0\mathrm{,rel}}(\hat{\vb{r}}_{ik})+U(\hat{\vb{r}}_{ik})$
with the non-interacting reduced Hamiltonian
$\ham_{0\mathrm{,rel}}(\hat{\vb{r}}_{ik})=\hbar^{2}\nabla^{2}_{ik}/2m^{*}$ and
reduced mass $m^{*}=m/2$. Splitting the Hamiltonian in this way, we
approximate the density matrix using the _pair product approximation_
Ceperley1995; Yan2015; Barker1979 given by

$\displaystyle\rho\left(\mathbf{R},\mathbf{R}^{\prime};\tau\right)\approx\left(\prod_{i=1}^{N}\rho_{1}\left(\mathbf{r}_{i},\mathbf{r}_{i}^{\prime},\tau\right)\right)\left(\prod^{N-1}_{i=i}\prod_{k=i+1}^{N}\frac{\rho^{\mathrm{rel}}\left(\mathbf{r}_{ik},\mathbf{r}_{ik}^{\prime};\tau\right)}{\rho_{0}^{\mathrm{rel}}\left(\mathbf{r}_{ik},\mathbf{r}_{ik}^{\prime};\tau\right)}\right)+\mathcal{O}(\tau^{3}),$
(24)

where $\rho_{1}$ is the _single-particle density matrix_ of the Hamiltonian
$\ham_{1}$, $\rho^{\mathrm{rel}}$ is the _two-body density matrix_ of the
interacting system $\ham^{\mathrm{rel}}$, and $\rho_{0}^{\mathrm{rel}}$ the
two-body density matrix of the non-interacting part. The form of the relative
contribution comes from the use of the Feynman-Kac formula FeyKac1949. The
exact derivation of the Feynman-Kac formula and pair product approximation can
be found in militzer2000path. In the following two subsections, I would like
to dissect and compute each part of the above expression.

##### Non-interacting propagator

The non-interacting single-particle density matrix $\rho_{1}$ is handled by
applying the primitive approximation as done in section 1. However, instead of
using BCH only ones, you can apply it recursively to generate higher-order
corrections, e.g.

$\displaystyle
e^{-\tau(\hat{A}+\hat{B})}=e^{-\frac{1}{2}\tau\hat{B}}e^{-\tau\hat{A}}e^{-\frac{1}{2}\tau\hat{B}}+O\left(\tau^{3}\right).$
(25)

This scheme of computing higher order is called _Trotter-Suzuki factorization_
suzuki_generalized_1976. In practice higher orders of Trotter-Suzuki
factorization are used, nevertheless, the second order is sufficient to
illustrate our purposes. Applying equation 25 to the non-interacting single
particle density matrix $\rho_{1}$, i.e. setting $\hat{A}=\ham_{0}$ and
$\hat{B}=V$, we can write

$\displaystyle\rho_{1}\left(\mathbf{R},\mathbf{R}^{\prime};\tau\right)=e^{-\frac{\tau}{2}\sum^{N}_{i=1}\left(V\left(\mathbf{r}_{i}\right)+V\left(\mathbf{r}_{i}^{\prime}\right)\right)}\mel{\mathbf{R}}{e^{-\tau\hat{H}_{0}}}{\mathbf{R}^{\prime}}.$
(26)

The only remaining unknown is the free particle propagator
$\rho_{0}\left(\mathbf{R},\mathbf{R}^{\prime};\tau\right)=\mel{\mathbf{R}}{e^{-\tau\hat{H}_{0}}}{\mathbf{R}^{\prime}}$.
Nevertheless, the free part of the Hamiltonian is $\hat{H}_{0}$ not diagonal
in the spatial basis. This can be resolved by decomposing the position
eigenstates in terms of the free particle plane waves or momentum eigenstates.
More concretely, we write

$\displaystyle\rho_{0}\left(\mathbf{R},\mathbf{R}^{\prime};\tau\right)$
$\displaystyle=\prod^{N}_{i=1}\int\frac{\dd{\vb{k}_{i}}}{(2\pi)^{d}}\int\frac{\dd{\vb{k}_{i}^{\prime}}}{(2\pi)^{d}}\mathrm{e}^{-i\vb{k}_{i}\cdot\vb{r}_{i}}\mathrm{e}^{i\vb{k}_{i}^{\prime}\cdot\vb{r}_{i}^{\prime}}\mel{\vb{k}_{1},\cdots,\vb{k}_{N}}{\mathrm{e}^{-\tau\lambda\sum_{j=1}^{N}\hat{\nabla}_{j}^{2}}}{\vb{k}_{1}^{\prime},\cdots,\vb{k}_{N}^{\prime}},$
(27)
$\displaystyle=\prod^{N}_{i=1}\int\frac{\dd{\vb{k}_{i}}}{(2\pi)^{d}}\int\frac{\dd{\vb{k}_{i}^{\prime}}}{(2\pi)^{d}}\mathrm{e}^{-\lambda\tau\left|\vb{k}_{i}^{\prime}\right|^{2}-i\vb{k}_{i}\cdot\vb{r}_{i}+i\vb{k}_{i}^{\prime}\cdot\vb{r}_{i}^{\prime}}\braket{\vb{k}_{1},\cdots,\vb{k}_{N}}{\vb{k}_{1}^{\prime},\cdots,\vb{k}_{N}^{\prime}}$
$\displaystyle=\prod^{N}_{i=1}\int\frac{d\vb{k}_{i}}{(2\pi)^{d}}\mathrm{e}^{-\lambda\tau\left|\vb{k}_{i}\right|^{2}+i\vb{k}_{i}\cdot\left(\vb{r}_{i}^{\prime}-\vb{r}_{i}\right)},$
(28)
$\displaystyle=(4\pi\lambda\tau)^{-Nd/2}\mathrm{e}^{-\frac{1}{4\lambda\tau}\sum_{i=1}^{N}\left(\vb{r}_{i}-\vb{r}_{i}^{\prime}\right)^{2}}.$
(29)

The result is the well-known expression for the free one particle propagator
in $d$ dimensions,

$\displaystyle\rho_{0}\left(\mathbf{r}_{i},\mathbf{r}_{i}^{\prime},\tau\right)=\frac{1}{\lambda_{\tau}^{d}}\exp\left(-\frac{\pi}{\lambda_{\tau}^{2}}\sum_{i=1}^{N}\left(\mathbf{r}_{i}-\mathbf{r}^{\prime}_{i}\right)^{2}\right),$
(30)

where $\lambda_{\tau}=\sqrt{4\pi\lambda\tau}$ is the de Broglie wavelength at
the temperature set by imaginary time step $\tau$. Substituting this in the
non-interacting single particle density matrix (26), we find

$\displaystyle\rho_{1}\left(\mathbf{R},\mathbf{R}^{\prime};\tau\right)=\frac{1}{\lambda_{\tau}^{dN}}\exp\left(-\frac{\pi}{\lambda_{\tau}^{2}}\sum_{i=1}^{N}\left(\mathbf{r}_{i}-\mathbf{r}^{\prime}_{i}\right)^{2}-\frac{\tau}{2}\sum^{N}_{i=1}\left(V\left(\mathbf{r}_{i}^{\prime}\right)+V\left(\mathbf{r}_{i}\right)\right)\right).$
(31)

Examining the result, we can already infer two findings. On the one hand,
smooth and straight paths generally have a higher probability of being
generated. On the other hand, paths with high potential energy have a low
probability of occurring.

There is a small detail I left out until now. Often, one wants to work with
periodic boundary conditions. This is ideal if one wants to study the bulk
properties of systems featuring infinite spatial extension. In addition, some
physical observables, such as superfluid fractions, are easier to measure if
one applies periodic boundary conditions to the simulation box. However, in
enforcing this, the size of the box becomes a key element in the numerical
approach, necessitating cautious extrapolation of all results. Consequently,
to write the momentum eigenstate in equation 2 over an integral, one must
assume that the wavelength of one imaginary time step $\tau$ is much smaller
than the size of the box with volume $V=L^{d}$, i.e., Ceperley1995;
militzer2000path; Spada2022

$\displaystyle\lambda_{\beta}\ll L.$ (32)

Not respecting this condition will lead to biased results. Recently, a new
formulation of WAPIMC has been published by Spada2022 which does not have such
a limitation111111Unfortunately, this was published way too late into my
thesis, but shows how relevant WAPIMC still is..

##### Interaction propagator

The interaction part is completely described by the second factor in the
expression (24) where we have that the numerator and denominator,
respectively, are given by

$\displaystyle\rho^{\mathrm{rel}}\left(\mathbf{r}_{ik},\mathbf{r}_{ik}^{\prime};\tau\right)=\mel{\mathbf{r}_{ik}^{\prime}}{\mathrm{e}^{-\tau\hat{H}^{\mathrm{rel}}}}{\mathbf{r}_{ik}}\qq{and}\rho^{\mathrm{rel}}_{0}\left(\mathbf{r}_{ik},\mathbf{r}_{ik}^{\prime};\tau\right)=\mel{\mathbf{r}_{ik}^{\prime}}{\mathrm{e}^{-\tau\hat{H}^{\mathrm{rel}}_{0}}}{\mathbf{r}_{ik}}.$
(33)

The latter, or the relative free density matrix $\rho^{\mathrm{rel}}_{0}$,
will be the same as expression (30) but with the replacement $m\to m^{*}=m/2$.
The former is often not that easy to find an expression for. The relative two-
body density matrix for a specific interaction potential can be computed by
solving the radial Schrödinger equation (see appendix
LABEL:app:_rel_int_prop). Nevertheless, often one has to perform a partial
wave decomposition, which leads to a non-closed analytical expression and
thereby must be numerically evaluated121212There exist other numerical methods
as well, such as the matrix squaring technique Ceperley1995.. However, in some
cases, the relative interaction propagator can be found in a closed analytical
form, e.g., for the one-dimensional delta or contact interaction potential
$U=g\delta\left(\mathbf{r}_{ij}\right)$ with $g$ an interaction coupling
constant Yan2015 computed

$\displaystyle\rho^{\mathrm{rel}}\left(\mathbf{r}_{ik},\mathbf{r}_{ik}^{\prime};\tau\right)=1-\exp\left(-\frac{m^{*}\left(\mathbf{r}_{ik}\mathbf{r}_{ik}^{\prime}+\left|\mathbf{r}_{ik}\mathbf{r}_{ik}^{\prime}\right|\right)}{\tau\hbar^{2}}\right)\times\sqrt{\frac{\pi
m^{*}\tau}{2}}\frac{g}{\hbar}\operatorname{erfc}(u)\exp(u),$ (34)

with
$u=m^{*}\left(r_{ik}+r_{ik}^{\prime}+g\tau\right)/\sqrt{2m^{*}\tau\hbar^{2}}$
and $\erfc(.)$ is the complementary error function. Another, noticeable
example is by CaoBerne1992 which found a three-dimensional closed analytical
approximation for the hard-sphere potential

$\displaystyle\frac{\rho^{\mathrm{rel}}\left(\mathbf{r}_{ik},\mathbf{r}_{ik}^{\prime};\tau\right)}{\rho^{\mathrm{rel}}_{0}\left(\mathbf{r}_{ik},\mathbf{r}_{ik}^{\prime};\tau\right)}=1-\frac{a\left(r_{ik}+r_{ik}^{\prime}\right)-a^{2}}{r_{ik}r_{ik}^{\prime}}e^{-\left[r_{ik}r_{ik}^{\prime}+a^{2}-a\left(r_{ik}+r_{ik}^{\prime}\right)\right](1+\cos\theta)m^{*}/\left(2h^{2}\tau\right)},$
(35)

with $a$ the cut-off (equivalent to the scattering length) of the hard-sphere
potential and $\theta$ the angle between $\vb{r}_{ik}$ and
$\vb{r}^{\prime}_{ik}$. To our knowledge, no two-dimensional closed analytical
form exist for any non-trivial interaction potential. In section
LABEL:sec:_derivation_int_prop we shall derive non-closed approximation which
can be numerically evaluated.

###  3 Monte Carlo Updates

In this section, we will construct the updates or moves to generate new
entries in the Markov chain. This is done by modifying the worldlines in the
configuration. One important thing to keep in mind is that we have to
guarantee the ergodicity of the updates, i.e., the whole configuration space
must be accessible starting from any configuration. In addition, any update
must respect the detailed balance equation (9).

To maintain the high acceptance rate and, thereby, make the move more
efficient, we should suggest new configurations in a way that the ratio in the
acceptance probability $\mathcal{A}$, i.e., equation 10, is around unity. One
way of doing this is to employ local updates, which means we only change a
small portion of the configuration such that the improbable configurations are
avoided. It is important to note that the efficiency of an update does not
affect the converged distribution, i.e., the obtained physics. It only affects
the speed of the convergence.

Let us start from a configuration $\config$. A simple but naive example for an
update would be moving a bead of worldline $i$ at time slice $j$ inside a
$d$-dimensional Ball $\mathcal{B}(\vb{r}_{i,j},\Delta)$ with radius $\Delta$
to a position $\vb{r}^{\prime}_{i,j}$.131313Note that the interval must
symmetric, as otherwise the detail balance is violated. The exact distance can
be uniformly generated from the cube. However, as discussed in section 2, in
the case one works with periodic boundary conditions, there can exist
ambiguity for the periodic image of the new position if $\Delta>L/2$ with $L$
being the length of the sides of the simulation box. Hence, one must ensure
$\Delta<L/2$.

The proposal probability
$\mathcal{P}\left(\config\rightarrow\config^{\prime}\right)$ to start from
configuration $\config$ and generate this new configuration $\config^{\prime}$
is than equal to starting form $\config^{\prime}$ and going to $\config$.
Hence, the acceptance probability defined in expression (10) is given by

$\displaystyle
A(\mathbf{X}\rightarrow\tilde{\mathbf{X}})=\min\left(1,e^{\tau\Delta
S(\mathbf{X}\rightarrow\tilde{\mathbf{X}})}\frac{\rho_{0}\left(\mathbf{r}_{i,j-1},\mathbf{r}^{\prime}_{i,j};\tau\right)\rho_{0}\left(\mathbf{r}^{\prime}_{i,j},\mathbf{r}_{i,j+1};\tau\right)}{\rho_{0}\left(\mathbf{r}_{i,j-1},\mathbf{r}_{i,j};\tau\right)\rho_{0}\left(\mathbf{r}_{i,j},\mathbf{r}_{i,j+1};\tau\right)}\right).$
(36)

Because of the $\beta$-periodicity we assume, without sacrificing generality,
throughout this thesis that whenever we write the $j+m$ time slice with
$m\in\mathbb{Z}$, the enumeration happens in modular $J$, i.e, $j+m\equiv
j+m\mod J$. Furthermore, notice, that the expression (10) is greatly reduced,
as only the affected time slices of the update contribute to the probability.
Nevertheless, the average acceptance probability will be low as for large
distances the kinetic free density matrices will vanish when the distance
between $\mathbf{r}_{i,j-1}$ or $\mathbf{r}_{i,j+1}$ and
$\mathbf{r}^{\prime}_{i,j}$ will become too large.

To make the updates more efficient, we make use of the fact that the proposal
probability $\mathcal{P}$ can be chosen freely. We can construct an update
such that the proposal probability $\mathcal{P}$ cancels factors in the
acceptance probability $\mathcal{A}$ to achieve unity for both the update and
the inverse. In other words, we perform _importance sampling_ , to generate a
configuration with a high acceptance probability.

For example, to expand on the above update, we could generate the new position
$r^{\prime}_{i,j}$ out of the kinetic propagators, such that the proposal
probability is given by
$\mathcal{P}\left(\config\rightarrow\config^{\prime}\right)=\tfrac{1}{JN}\rho_{0}\left(\mathbf{r}_{i,j-1},\mathbf{r}^{\prime}_{i,j};\tau\right)\rho_{0}\left(\mathbf{r}^{\prime}_{i,j},\mathbf{r}_{i,j+1};\tau\right)$
and the acceptance probability becomes141414The reverse proposal probability
is given by
$\mathcal{P}\left(\config^{\prime}\rightarrow\config\right)=\tfrac{1}{JN}\rho_{0}\left(\mathbf{r}_{i,j-1},\mathbf{r}_{i,j};\tau\right)\rho_{0}\left(\mathbf{r}_{i,j},\mathbf{r}_{i,j+1};\tau\right)$.

$\displaystyle A(\config\to\config^{\prime})=\min\left(1,e^{\tau\Delta
S(\config\to\config^{\prime})}\right).$ (37)

The $\tfrac{1}{JN}$ comes from selecting a bead to move. This is possible
because the kinetic propagator can be seen as a Gaussian distribution. There
are only a couple of probability distributions where one can sample directly,
and the Gaussian is one of them. The remaining factor in the above expression
cannot be removed, as one cannot simply sample a new configuration out of it.

There is one last problem to consider. Single time slice moves are extremely
insufficient. Not because the acceptance probability is small, but because the
newly subsequent generated configurations are so similar that the
configuration space is explored too slowly. Ceperley1995 found that the amount
of computation time required to modify the overall form of a single path
scales as $J^{3}$. This will be addressed in the next subsections, where we
formulate the official updates used in our WAPIMC code.

####  1 Canonical updates

##### Center-of-mass update

Instead of moving individual beads, one could just move a whole worldline by a
distance $\vb{r}^{\prime}$. Now, remember that the worldlines are permuted.
Hence, one can not simply move one worldline, as the link connecting the last
and first bead will become highly improbable, and thereby, the update gets
easily rejected in the Metropolis-Hastings question. To solve this one has to
move to the whole permutation cycle or polymer. The move is called the
_centre-of-mass update_ as only the centre-of-mass of the polymer is moved —
the relative position of the bead on the polymer remains the same. This has
the consequence that there will be no kinetic component in the acceptance
probability. The move is very efficient if the average polymer size is less
than the polymer distance or the temperature is higher than the degeneracy
temperature. This is because there will be less overlap between the polymers.
At low temperatures, it can be more efficient to only move one particle of the
polymer.

Figure 2: Illustration of the Center-of-mass update. A permutation cycle of
worldlines or polymer gets translated by a vector $\vb{r}^{\prime}$. Notice
the links between two consecutive beads are not modified. The unchanged
worldlines are portrayed in black, the old beads and links are dashed, and the
modified moved closed worldline $i$ is in solid grey.

More concretely, the centre-of-mass update goes as follows:

  1. [i)]

  2. 1.

Choose a single worldline $i$ uniformly from the last configuration in the
Markov chain.

  3. 2.

Sample the polymer displacement $\vb{r}^{\prime}$ uniformly from the symmetric
interval $\left[-\Delta,\Delta\right]^{d}$. $\Delta$ is a free parameter that
can be optimized to have a desirable acceptance ratio of the centre-of-mass
update during the simulation.

  4. 3.

Move the polymer by adding $\mathbf{r}^{\prime}$ to the position of all the
beads and compute the weights of the old and new configurations.

  5. 4.

Determine if the new configuration is accepted by the Metropolis-Hastings
question with the acceptance probability

$A_{\text{Move
}}(\config\rightarrow\config^{\prime})=\min\left(1,e^{\tau\Delta
S(\config\rightarrow\config^{\prime})}\right).$

##### Reshape update

Let us go back to modifying individual beads. Instead of moving one individual
bead of a polymer by sampling the new position out of the kinetic propagator
$\rho_{0}$, we could modify $M-2$ beads of the polymer starting from a time
slice $j_{0}$ and ending at $j_{M}$, with $M$ smaller than the total length of
the selected polymer. This procedure of generating the path
$\left(\mathbf{r}_{j_{0}},\mathbf{r}_{j_{0}+1}^{\prime},\ldots,\mathbf{r}_{j_{0}+M-1}^{\prime},\mathbf{r}_{j_{M}}\right)$
is called generating a _Brownian bridge_ krauth2006 between $j_{0}$ and
$j_{M}$. The proposal probability of the direct and reverse update is given by

$\displaystyle\mathcal{P}\left(\config\rightarrow\config^{\prime}\right)$
$\displaystyle=\frac{1}{JN}\rho_{0}\left(\mathbf{r}_{\mathrm{p}(j_{0}),j_{0}},\mathbf{r}^{\prime}_{\mathrm{p}(j_{0}+1),j_{0}+1};\tau\right)\cdots\rho_{0}\left(\mathbf{r}^{\prime}_{\mathrm{p}(j_{M}-1),j_{M}-1},\mathbf{r}_{\mathrm{p}(j_{M}),j_{M}};\tau\right),$
(38) $\displaystyle\mathcal{P}\left(\config^{\prime}\rightarrow\config\right)$
$\displaystyle=\frac{1}{JN}\rho_{0}\left(\mathbf{r}_{\mathrm{p}(j_{0}),j_{0}},\mathbf{r}_{\mathrm{p}(j_{0}+1),j_{0}+1};\tau\right)\cdots\rho_{0}\left(\mathbf{r}_{\mathrm{p}(j_{M}-1),j_{M}-1},\mathbf{r}_{\mathrm{p}(j_{M}),j_{M}};\tau\right),$
(39)

with $\mathrm{p}(j)$ the worldline index of the particle in the polymer at
time slice $j$. The resulting acceptance probability is then

$\displaystyle A(\config\to\config^{\prime})=\min\left(1,e^{\tau\Delta
S(\config\to\config^{\prime})}\right).$ (40)

In figure 3, a visualization of such a procedure is shown in the left
worldline. Notice the Gaussian distributions at every newly sampled bead
centred around the previous bead. This is the Gaussian where the bead is
sampled from. However, it only depends on the previous bead and is thus
independent of the last bead $\mathbf{r}_{j_{M}}$ of the new path. This has
the result that the last generated bead can be far from $\mathbf{r}_{j_{M}}$,
which leads to a very small acceptance probability.

Figure 3: Illustration of the different reshape schemes. The left path is
sampled using the proposal naive probability whereas on the right the path is
sampled taking the endpoint into account. The Gaussian distributions are shown
for the first three beads of every Brownian bridge. The new paths are
portrayed in solid grey.

A solution would be to generate the Brownian bridge with the Levy construction
krauth2006, which takes the last bead $\mathbf{r}_{j_{M}}$ of the path into
account. This adds up to complementing each factor in the proposal probability
(38) with a kinetic propagator to the last bead $\mathbf{r}_{j_{M}}$, which
yields

$\mathcal{P}\left(\config\rightarrow\config^{\prime}\right)=\frac{1}{JN}\frac{\rho_{0}\left(\mathbf{r}_{j_{0}},\mathbf{r}^{\prime}_{j_{0}+1},\tau\right)\rho_{0}\left(\mathbf{r}^{\prime}_{j_{0}+1},\mathbf{r}_{j_{M}},(M-1)\tau\right)}{\rho_{0}\left(\mathbf{r}_{j_{0}},\mathbf{r}_{j_{M}},M\tau\right)}\times\\\
\hskip
28.45274pt\frac{\rho_{0}\left(\mathbf{r}^{\prime}_{j_{0}+1},\mathbf{r}^{\prime}_{j_{0}+2},\tau\right)\rho_{0}\left(\mathbf{r}^{\prime}_{j_{0}+2},\mathbf{r}_{j_{M}},(M-2)\tau\right)}{\rho_{0}\left(\mathbf{r}^{\prime}_{j_{0}+1},\mathbf{r}_{j_{M}},(M-1)\tau\right)}\times\cdots\\\
\cdots\times\frac{\rho_{0}\left(\mathbf{r}^{\prime}_{j_{M}-2},\mathbf{r}^{\prime}_{j_{M}-1},\tau\right)\rho_{0}\left(\mathbf{r}^{\prime}_{j_{M}-1},\mathbf{r}_{j_{M}},\tau\right)}{\rho_{0}\left(\mathbf{r}^{\prime}_{j_{M}-2},\mathbf{r}_{j_{M}},2\tau\right)}.$
(44)

Here, we neglected the index notation for the worldlines to make the
expression more readable. We can simplify the expression by cancelling some
factors such that the direct and reverse proposal probabilities give

$\displaystyle\mathcal{P}_{\mathrm{Reshape}}\left(\config\rightarrow\config^{\prime}\right)$
$\displaystyle=\frac{1}{JN}\frac{\rho_{0}\left(\mathbf{r}_{\mathrm{p}(j_{0}),j_{0}},\mathbf{r}^{\prime}_{\mathrm{p}(j_{0}+1),j_{0}+1},\tau\right)\cdots\rho_{0}\left(\mathbf{r}^{\prime}_{\mathrm{p}(j_{M}-1),j_{M}-1},\mathbf{r}^{\prime}_{\mathrm{p}(j_{M}),j_{M}},\tau\right)}{\rho_{0}\left(\mathbf{r}_{\mathrm{p}(j_{0}),j_{0}},\mathbf{r}_{\mathrm{p}(j_{M}),j_{M}},M\tau\right)},$
(45)
$\displaystyle\mathcal{P}_{\mathrm{Reshape}}\left(\config^{\prime}\rightarrow\config\right)$
$\displaystyle=\frac{1}{JN}\frac{\rho_{0}\left(\mathbf{r}_{\mathrm{p}(j_{0}),j_{0}},\mathbf{r}_{\mathrm{p}(j_{0}+1),j_{0}+1},\tau\right)\cdots\rho_{0}\left(\mathbf{r}_{\mathrm{p}(j_{M}-1),j_{M}-1},\mathbf{r}_{\mathrm{p}(j_{M}),j_{M}},\tau\right)}{\rho_{0}\left(\mathbf{r}_{\mathrm{p}(j_{0}),j_{0}},\mathbf{r}_{\mathrm{p}(j_{M}),j_{M}},M\tau\right)},$
(46)

and the acceptance probability of equation 40 remains the same. We call the
move corresponding to the above proposal probabilities the _reshape update_.
In practice, a position $\vb{r}_{k}$ of a single bead in the path with
$k\in[1,M]$ is sampled out of the Gaussian probability distribution per
dimension

$\displaystyle
P\left(r_{k}\right)=\frac{1}{\sqrt{2\pi}\sigma_{k}}\exp\left(-\frac{\left({r}_{k}-{\xi}_{k}\right)^{2}}{2{\sigma}_{k}^{2}}\right),$
(47)

with the variance and mean, respectively, given by

$\displaystyle{\sigma}_{k}=\sqrt{\frac{\alpha_{k}}{2\pi}}\lambda_{\tau}=\sqrt{2\lambda\alpha_{k}\tau}\qq{and}{\xi}_{k}=\alpha_{k}{r}_{k-1}+\left(1-\alpha_{k}\right){r}_{M+1},\qq{where}\alpha_{k}=\frac{M+1-k}{M+2-k}.$
(48)

In figure 3, the difference in the sampling is shown in the right worldline.
Now, the Gaussian where the beads are sampled from is not centred around the
previous bead any more. Instead, the intersection of the line between the
previous bead and the last bead of the path with the horizontal line of the
time slice serves as the mean of the Gaussian. The advantage being that the
last bead of the path is taken into account, and thus large gaps, as in figure
3 on the right, do not occur.

In conclusion, the reshape update is formulated as follows:

  1. [i)]

  2. 1.

Choose a single worldline uniformly from the last configuration in the Markov
chain.

  3. 2.

Select the time slice $j_{0}\in[1,J]$ and number of beads to be reshaped
$M\in[1,\bar{M}]$ uniformly. $\bar{M}$ is a free parameter to be optimized
during the WAPIMC simulation.

  4. 3.

Construct the new path from time slice $j_{0}$ until $j_{M}=j_{0}+M$ with the
Levy construction.

  5. 4.

Determine if the new configuration is accepted by the Metropolis-Hastings
question with the acceptance probability given by equation 40.

##### Permutation update

Until now, no update has been introduced which can sample configurations where
permutations between worldlines are present. However, these permutations
represent the quantum statistics of the system, which is a crucial ingredient
of the physics of ensembles of indistinguishable particles. It underlies
phenomena such as Bose-Einstein condensation and superfluidity — both
conditions of interest for our purposes.

With the reshape update in mind, one can easily think of a similar update that
can permute the worldlines. Instead of creating a Levy construction between
two beads of the same worldlines, one instead creates a bridge between beads
of two different worldlines. A simple example and comparison with the reshape
update is shown in figure 4. The same proposal probabilities analysis as in
the reshape update can be worked out, which results in the same acceptance
probability as in equation 40.

Naturally, the update is likely to be rejected if the to-be connected beads of
the corresponding two worldlines are far apart. Hence, one can again do
importance sampling. Suppose we have selected a worldline $i_{1}$. We can
sample a worldline $i_{2}$ in the neighbourhood of $i_{1}$ with the
probability

$\displaystyle
p\left(i_{2}\right)=\frac{\rho_{0}\left(\mathbf{r}_{i_{1},j_{0}},\mathbf{r}_{i_{2},j_{M}};M\tau\right)\rho_{0}\left(\mathbf{r}_{i_{2},j_{0}},\mathbf{r}_{i_{1},j_{M}};M\tau\right)}{\sum_{i=1}^{N}\rho_{0}\left(\mathbf{r}_{i_{1},j_{0}},\mathbf{r}_{i,j_{M}};M\tau\right)\rho_{0}\left(\mathbf{r}_{i,j_{0}},\mathbf{r}_{i_{1},j_{M}};M\tau\right)}.$
(49)

This way we guarantee that the distance
$\abs{\mathbf{r}_{i_{1},j_{0}}-\mathbf{r}_{i_{2},j_{M}}}$ is of the order
$\sqrt{M\tau}$, consequently, the permutation update is more probable to be
accepted, and the sampling process can proceed more swiftly across the
configuration space.

In conclusion, the permutation update is formulated as follows:

  1. [i)]

  2. 1.

Choose a single worldline uniformly from the last configuration in the Markov
chain.

  3. 2.

Select the time slice $j_{0}\in[1,J]$ and number of beads to be reshaped
$M\in[1,\bar{M}]$ uniformly. $\bar{M}$ is a free parameter to be optimized
during the WAPIMC simulation.

  4. 3.

Sample a second wordline $i_{2}$ from the probability distribution (49). If
$i_{1}=i_{2}$ reject the update.

  5. 4.

Construct the new paths from the bead $(i_{1},j_{0})$ to
$(i_{2},j_{M}=j_{0}+M)$ and bead $(i_{2},j_{0})$ to $(i_{1},j_{M}=j_{0}+M)$
with the Levy construction.

  6. 5.

Determine if the new configuration is accepted by the Metropolis-Hastings
question with the acceptance probability given by equation 40.

There is still one major problem with this type of permutation sampling. In
the presence of repulsive and/or hard-core potential, the permutation update
is bound to become inefficient, i.e., the proposed update has a high
likelihood of being rejected. This is because when there are a lot of
particles, it becomes impossible to sample a path in an unbiased way where no
two particles become close to each other. As a result, it is hard to generate
long exchange cycles, which are the configurations related to the superfluid
fraction, and thereby, lead to large errors151515It will lead to large
autocorrelation times, which is translated into an error (see 5).. This is an
ergodicity problem as the whole configuration space can not be sampled in a
reasonable computation time. This worsens exponentially in the function of the
number of particles $N$. In general, such a problem is called a _critical
slowing down phenomenon_.

Figure 4: Illustration of the reshape red(a) and permutation red(b) update
with the levy construction. The unchanged beads are portrayed in black, the
old beads and links are dashed, and the new path or Brownian bridge is
depicted in solid grey.

####  2 Worm algorithm

[remember picture, overlay] [opacity=1,inner sep=0pt] at ($(0,0)-(0cm,16cm)$);
To deal with this critical slowing down phenomenon of generating permuted
configurations, the Worms Algorithm was created, initially in the discrete
spatial variation of PIMC prokofev1998sec and afterwards in the continuous
version Boninsegn2006PRL; Boninsegn2006. Figure 5 depicts the primary concept.
Certain sections of the configuration space are difficult to reach with the
help of local updates described above, so, why not extend the configuration
space by allowing _open worldlines_? In that way, one could take a detour in
this so-called _G-sector_ 161616We call this subset of the grand-canonical
configuration space the G-sector because the configuration is associated with
sampling the Matsubara Green’s function., i.e., the set of all configurations
in the ensemble that contain open worldlines or _worms_ , to generate
configurations in the complement or Z-sector. In figure 5, an example of a
configuration in the G-sector is shown. For historical reasons, the tail of
the worm, i.e., the first bead, is called _Masha_ $\mathcal{M}$ and the head
or last bead is called _Ira_ $\mathcal{I}$.

[b].45 [width=0.85]Chapter_4/Figures/extendingconfig.pdf (a)[b].45
[width=0.90]Chapter_4/Figures/worm.pdf (b)

Figure 5: Figures adjoining the Worm Algorithm. fig: extendingconfig Sketch
representing a walk in the configuration space from an initial configuration
$\config$ to a final configuration $\config^{\prime}$. The walker takes a
detour in the G-sector to reach a point in the Z-sector more easily. fig:
Gsector config A configuration in the G-sector with one closed worldline and a
worm. The head (Ira) and tail (Masha) are labelled by $\ira$ and $\masha$.

The worm, or open worldline, of the G-sector configurations break the critical
slowing down phenomena by serving as a mediator and assisting the worldlines
to generated permutations by using the Monte Carlo update called _swap_ (see
section 2). Furthermore, extending the configuration space allows PIMC to
sample the grand-canonical ensemble. This is best explained with the help of
figure 6. One introduces Monte Carlo updates, called _advance_ and _recede_ ,
so the worm can, respectively, grow and shrink. By virtue of the
$\beta$-periodicity of the worldlines, the particle number or closed wordlines
can change if one back closed the worm.

Figure 6: Four worldline configurations demonstrating how the number of closed
worldlines $N$ can increase with the help of the Worm Algorithm updates. First
the closed worldine is opened to a worm (a $\to$ b). Next, the head of the
worm is advanced (b $\to$ c). Lastly, the worm is closed (c $\to$ d). For each
state the number of closed worldlines $N$ is given.

Using the Worm algorithm has other advantages other than solving ergodicity
problem for permutated configurations. It provides PIMC with additional
physical quantities (see section 4). In addition to accurately determining the
superfluid fraction $f$, one can determine the compressibility $\kappa$, as it
relates to the fluctuation of the number of particles. Furthermore, as one
also generates configuration in the G-sector, one can measure observables
related to the _Matsubara Green’s function_ , such as correlation functions
and the momentum distribution. Lastly, also dynamical properties can be
extracted, as the Matsubara Green function is connected to the single-particle
excitation spectrum through a Laplace transform dornheim2014;
FilinovBonitz2012.

##### Extending the partition function

In this section, the formalism of the Worms Algorithm is explained. It is
represented in the grand canonical ensemble. Nevertheless, the worm algorithm
in the canonical ensemble is done analogously by taking the limit $\mu\to 0$.

First, let us introduce some jargon. As alluded to earlier, the configuration
space of the grand-canonical ensemble is split into two complemented parts or
sectors, the G-sector and Z-sector, which refer to the configuration
containing a worm or _off-diagonal_ configuration and the closed configuration
or _diagonal_ configurations, respectively.

To add the Worms Algorithm to PIMC we must rewrite the partition function in
terms of the grand-canonical Hamiltonian $\hat{K}=\ham-\mu\hat{N}$ where
$\ham$ is the Hamiltonian of the canonical ensemble, $\hat{N}$ is the particle
number operator, and $\mu$ the chemical potential. A direct consequence will
be that the particle number $N$ can fluctuate and instead the chemical
potential $\mu$ will be the third fixed-parameter besides the volume and
temperature. Given the definition of the partition function (1), we can write

$\displaystyle
Z_{\mathrm{GC}}=\sum_{N=0}^{+\infty}\frac{1}{N!}\sum_{\sigma\in\Pi}\int\dd{\vb{R}}\mel{\sigma\cdot\mathbf{R}}{e^{-\beta\hat{K}}}{\mathbf{R}},$
(50)

where the trace was expanded in the spatial position basis and all the
possible particle numbers. Again, as in section 1, splitting $\beta$ into $J$
imaginary time slices of size $\tau$ yields the weight of a closed
configuration with $N$ worldlines or particles

$\displaystyle\mathcal{W}(\config)=e^{\beta\mu
N}\mel{\sigma\cdot\mathbf{R}_{\mathbf{0}}}{e^{-\tau\hat{H}}}{\mathbf{R}_{\mathbf{J}-\mathbf{1}}}\times\ldots\times\mel{\mathbf{R}_{\mathbf{1}}}{e^{-\tau\hat{H}}}{\mathbf{R}_{\mathbf{0}}}.$
(51)

Notice the $1/N!$ factor disappears as there are $N!$ variations in labelling
the worldlines in the configuration.

As explained above, we do not only want to extend the configuration space to
the grand-canonical ensemble with closed worldlines but also add the
configuration containing open worldlines. The latter does not contribute to
the partition function (50), but rather to an additional term given by

$\displaystyle
Z^{\prime}=Z_{\mathrm{GC}}\sum_{j_{\masha}}\sum_{j_{\ira}}\int\dd{\vb{r}_{\masha}}\dd{\vb{r}_{\ira}}G(\vb{r}_{\masha},\vb{r}_{\ira};\tau(j_{\ira}-j_{\masha})),$
(52)

where the under scripts $\masha$ and $\ira$ stand for Masha and Ira,
respectively, and $G$ is the Matsubara Green’s function given by

$\displaystyle
G(\vb{r}_{j},\vb{r}_{k};\tau)=\expval{\hat{\mathcal{T}}\hat{\psi}(\vb{r}_{\ira},\tau)\hat{\psi}^{\dagger}(\vb{r}_{\masha},0)}.$
(53)

Here is $\hat{\mathcal{T}}$ the time-ordering operator, $\hat{\psi}$ and
$\hat{\psi}^{\dagger}$, respectively, the bosonic annihilation and creation
operators in the Matsubara representation. The latter are both operators
evolved in imaginary time in the Heisenberg representation
altland_simons_2010. The Green’s function corresponds to the worm as it
creates a particle at the bead $\masha$ and annihilates it at $\ira$. The
generalized partition function becomes
$Z_{\mathrm{tot}}=Z_{\mathrm{GC}}+CZ^{\prime}$ where $C\in\mathbb{R}^{+}$ is
an optimization parameter which remains constant during simulation. The
constant parameter determines the ratio of the diagonal and off-diagonal
configurations during the WAPIMC simulation. More concretely, by labelling the
number of diagonal configurations sampled in the Markov chain as
$N_{\mathrm{Z}}$ and off-diagonal configurations as $N_{\mathrm{G}}$, so that
the total number of configurations in the Markov chain
$N_{\mathrm{MC}}=N_{\mathrm{Z}}+N_{\mathrm{G}}$, we can write

$\displaystyle\frac{N_{\mathrm{Z}}}{N_{\mathrm{MC}}}=\frac{Z_{\mathrm{GC}}}{Z_{\mathrm{tot}}}\qq{and}\frac{N_{\mathrm{G}}}{N_{\mathrm{MC}}}=C\frac{Z^{\prime}}{Z_{\mathrm{tot}}}\qq{if}N_{\mathrm{MC}}\to\infty.$
(54)

Combining both expressions yields

$\displaystyle\frac{N_{\mathrm{G}}}{N_{\mathrm{Z}}}=C\frac{Z_{\mathrm{G}}}{Z}\qq{if}N_{\mathrm{MC}}\to\infty.$
(55)

It follows that $C$ is linearly proportional to the ratio
$N_{\mathrm{G}}/N_{\mathrm{Z}}$. This relation allows for equilibrated
transitions between Z and the G-sector.

##### Updates in the Worm Algorithm

The Worms Algorithm proposes an additional set of Monte Carlo Markov chain
updates that modify the worms of the configurations, thereby, changing the
configuration to the G- or Z-sector. Naturally, they must fulfil the detail
balanced (9) and ergodicity conditions. In the following, we present the set
updates for the grand-canonical ensemble originally constructed by
Boninsegn2006PRL. For each update, the proposal and acceptance probabilities
are discussed. Afterwards, the set of updates is adjusted to the sample in the
canonical ensemble.

  * [leftmargin=*]

  * •

Insert and close: This pair of complementary update “inserts” a worm into a
diagonal configuration or “removes” a worm from an off-diagonal configuration.
Hence, a transition from the Z-sector to the G-sector is made or vice versa.
Together, the updates respect the detailed-balance condition. An example can
be seen in figure 7. The insertion of a worm goes by the following steps:

    1. [i)]

    2. 1.

Choose the time slice $j_{0}=j_{\masha}\in[1,J]$ and spatial coordinate
$\vb{r}_{\masha}$ for the tail $\masha$ uniformly in the volume $V$ of the
system.

    3. 2.

Sample the number of beads to be created $M\in[1,\bar{M}]$ uniformly.
$\bar{M}$ is a free parameter to be optimized during the WAPIMC simulation.

    4. 3.

Sample the coordinate of the worm’s head on time slice
$j_{\ira}=j_{M}=j_{0}+M$ from the kinetic distribution
$\rho_{0}(j_{\masha},j_{\ira};M\tau)$.

    5. 4.

Construct the open worldline from the tail $\masha$ to the head $\ira$ with
the Levy construction by (47).

As a result, we get the proposal probability

$\displaystyle\mathcal{P}_{\mathrm{Insert}}(\config\to\config^{\prime})=\frac{\rho_{0}(\vb{r}_{\masha},\vb{r}_{\ira};M\tau)}{VJ\bar{M}}\prod^{j_{\ira}-1}_{j=j_{\masha}}\frac{\rho_{0}(\vb{r}_{j},\vb{r}_{j+1};\tau)}{\rho_{0}(\vb{r}_{\masha},\vb{r}_{\ira};M\tau)}=\frac{1}{VJ\bar{M}}\prod^{j_{\ira}-1}_{j=j_{\masha}}\rho_{0}(\vb{r}_{j},\vb{r}_{j+1};\tau),$
(56)

where denominator contribution from the levy construction
$\rho_{0}(j_{\masha},j_{\ira};M\tau)$ cancels because we sample the worm’s
head from the kinetic distribution $\rho_{0}(j_{\masha},j_{\ira};M\tau)$. The
removal goes as follows:

    1. [i)]

    2. 1.

Count the number of links $M$ the worm contains. If $M>\bar{M}$ reject the
update where $\bar{M}$ is the same optimization parameter as for insert.

    3. 2.

Remove the worldline.

The proposal probability of the remove update is given by
$\mathcal{P}_{\mathrm{Insert}}(\config^{\prime}\to\config)=1$. Hence, the
acceptance probabilities of the insert and remove update yield

$\displaystyle\mathcal{A}_{\mathrm{Insert}}(\config\to\config^{\prime})$
$\displaystyle=\min\left(1,CVJ\bar{M}e^{\tau M\mu}e^{\tau\Delta
S(\config\to\config^{\prime})}\right),$ (57)
$\displaystyle\mathcal{A}_{\mathrm{Remove}}(\config^{\prime}\to\config)$
$\displaystyle=\min\left(1,\frac{1}{CVJ\bar{M}}e^{-\tau M\mu}e^{-\tau\Delta
S(\config^{\prime}\to\config)}\right).$ (58)

It is important to again state that the optimization parameter $\bar{M}$ must
be the same for both insert and remove. Otherwise, the detailed-balance
condition is not satisfied.

Figure 7: Illustration of the insert red(a) and remove red(b) update. The
unchanged beads are portrayed in black, the old beads and links are dashed,
and the newly generated path or Brownian bridge in solid grey.

  * •

Open and Close: This pair of complementary update “opens” a closed worldline
in a diagonal configuration or “closes” the worm in an off-diagonal
configuration. Hence, a transition from the Z-sector to the G-sector is made
or vice versa. An example can be seen in figure 8. Together, the updates
respect the detailed-balance condition (9). The opening of a diagonal
worldline goes by the following steps:

    1. [i)]

    2. 1.

Choose uniformly a closed worldline $i$ from all $N$ worldlines and a time
slice $j_{0}=j_{\ira}\in[1,J]$ for the head $\ira$ of the new worm.

    3. 2.

Sample the number of links to be removed $M\in[1,\bar{M}]$ uniformly.
$\bar{M}$ is a free parameter to be optimized during the WAPIMC simulation.

Hence, the proposal probability is therefore given as

$\displaystyle\mathcal{P}_{\mathrm{open}}(\config\to\config^{\prime})=\frac{1}{JN\bar{M}}.$
(59)

The close update is implemented as follows:

    1. [i)]

    2. 1.

Count the number of beads $M$ that must be constructed to close the worm. If
$M>\bar{M}$ reject the update where $\bar{M}$ is the same optimization
parameter as for open.

    3. 2.

Close the open worldline from the head $\ira$ to the tail $\masha$ with the
Levy construction by (47).

which yields a proposal probability

$\displaystyle\mathcal{P}_{\mathrm{close}}(\config^{\prime}\to\config)=\prod^{j_{\masha}-1}_{j=j_{\ira}}\frac{\rho_{0}(j,j+1;\tau)}{\rho_{0}(j_{\masha},j_{\ira};M\tau)}.$
(60)

The above proposal probabilities lead to acceptance probabilities

$\displaystyle\mathcal{A}_{\text{Close }}(\config\to\config^{\prime})$
$\displaystyle=\min\left(1,\frac{\rho\left(\mathbf{r}_{\masha},\mathbf{r}_{\ira};M\tau\right)e^{\tilde{M}\tau\mu}e^{\tau\Delta
S(\config\to\config^{\prime})}}{CJN\bar{M}}\right),$ (61)
$\displaystyle\mathcal{A}_{\text{Open }}(\config^{\prime}\to\config)$
$\displaystyle=\min\left(1,\frac{CJN\bar{M}e^{-M\tau\mu}e^{-\tau\Delta
S(\config^{\prime}\to\config)}}{\rho\left(\mathbf{r}_{\masha},\mathbf{r}_{\ira};M\tau\right)}\right).$
(62)

Again, note that the optimisation parameter $\bar{M}$ must be the same for
both open and close when the update is proposed to advance the Markov chain.
Furthermore, $N$ is the number of closed worldlines if we open the
configuration. For example, if we want to propose to close a configuration in
the G-sector which only contains a worm with a bead smaller than $M$, the
reverse update would be to open the closed configuration, so $N=1$.

Figure 8: Illustration of the open red(a) and close red(b) update. The
unchanged beads are portrayed in black, the old beads and links are dashed,
and the newly generated path or Brownian bridge in solid grey.

  * •

Advance and recede: Of course, only modifying configurations to contain or
disclude a worm is not enough. We would also like to modify the beads of the
worm. Naturally, the canonical updates reshape and center-of-mass of section 1
can be used to adjust the open worldline. However, we would also like to grow
or shrink the worm. This is done with the advance and recede update,
respectively. An example is shown in figure 9. To advance the worm, we undergo
the steps:

    1. [i)]

    2. 1.

Sample the number of beads to be created $M\in[1,\bar{M}]$ uniformly.
$\bar{M}$ is a free parameter to be optimised during the WAPIMC simulation.

    3. 2.

Sample the coordinate of the new worm’s head on time slice
$j_{\ira}=j_{M}=j_{0}+M$ from the kinetic distribution
$\rho_{0}(\vb{r}_{0},\vb{r}_{\ira};M\tau)$ where $j_{0}$ and $\vb{r}_{0}$ are
the time slice and position of the old head, respectively.

    4. 3.

Construct the worldline from the old head to the new head with the Levy
construction by (47).

Conversely, to recede the worm we must:

    1. [i)]

    2. 1.

Sample the number of beads to be removed $M\in[1,\bar{M}]$ uniformly.
$\bar{M}$ is the same free parameter as for the advance update. If $M$ is
greater than the number of beads the worm contains, reject the update.

    3. 2.

Remove the beads of the worm on the time slices $j_{\ira}$, $\ldots$,
$j_{\ira-M}$.

The above procedures lead to the proposal probabilities

$\displaystyle\mathcal{P}_{\mathrm{Advance}}(\config\to\config^{\prime})$
$\displaystyle=\frac{\rho_{0}(\vb{r}_{0},\vb{r}_{\ira};M\tau)}{\bar{M}}\prod^{j_{\ira}-1}_{j=j_{\masha}}\frac{\rho_{0}(\vb{r}_{j},\vb{r}_{j+1};\tau)}{\rho_{0}(\vb{r}_{0},\vb{r}_{\ira};M\tau)},$
(63) $\displaystyle\mathcal{P}_{\mathrm{Recede}}(\config\to\config^{\prime})$
$\displaystyle=\frac{1}{\bar{M}}.$ (64)

As a result, the acceptance probabilities of the advance and recede update in
the Metropolis question yields

$\displaystyle\mathcal{A}_{\text{Advance}}(\config\to\config^{\prime})$
$\displaystyle=\min\left(1,e^{\tilde{M}\tau\mu}e^{\tau\Delta
S(\config\to\config^{\prime})}\right),$ (65)
$\displaystyle\mathcal{A}_{\text{Recede}}(\config^{\prime}\to\config)$
$\displaystyle=\min\left(1,e^{-M\tau\mu}e^{-\tau\Delta
S(\config^{\prime}\to\config)}\right).$ (66)

Once more, the optimisation parameter $\bar{M}$ must be the same for advance
and recede when one of them is proposed.

Figure 9: Illustration of the advance red(a) and recede red(b) update. The
unchanged beads are portrayed in black, the old beads and links are dashed,
and the newly generated path or Brownian bridge in solid grey.

  * •

Swap: Last but not least, we must address the critical slowing down phenomenon
where we introduce the Worm algorithm. That is the sampling of permutations
for the configuration in the Z-sector, which is done with the swap update. The
concept is best explained visually which can be found in figure 10. The steps
are as follows:

    1. [i)]

    2. 1.

The number of beads to be modified is $\bar{M}$. $\bar{M}$ is a free parameter
to be optimised during the WAPIMC simulation.

    3. 2.

Construct the list $\chi$ of all kinetic propagators
$\rho_{0}(\vb{r}_{\ira},\vb{r}_{k,j_{\bar{M}}},\bar{M}\tau)$ with
$k=1,\ldots,N$, $j_{\bar{M}}=j_{\ira}+\bar{M}$ and $N$ the number of closed
worldlines.

    4. 3.

Sample a bead $\alpha$ discretely by the weights $w(k)=\chi(k)/\Sigma$ where
$\Sigma=\sum^{N}_{k=1}\chi(k)$.

    5. 4.

Locate the bead $\delta$ on the time slice of the head of the worm $\ira$ by
backtracking on the polymer from the bead $\alpha$ with $\bar{M}$ beads. Label
the path as
$(\mathbf{r}_{j_{\ira}},\mathbf{r}_{j_{\ira}+1},\ldots,\mathbf{r}_{j_{\ira}+\bar{M}-1},\mathbf{r}_{j_{\ira}+\bar{M}})$
where $\mathbf{r}_{j_{\ira}}=\mathbf{r}_{\delta}$ and
$\mathbf{r}_{j_{\ira}+\bar{M}}=\mathbf{r}_{\alpha}$.

    6. 5.

Construct the list $\chi^{\prime}$ of all kinetic propagators
$\rho_{0}(\vb{r}_{\delta},\vb{r}_{k,j_{\bar{M}}},\bar{M}\tau)$ with
$k=1,\ldots,N$ and compute $\Sigma^{\prime}=\sum^{N}_{k=1}\chi^{\prime}(k)$.

    7. 6.

Remove all beads in between $\alpha$ and $\delta$ and generate a new path
$(\mathbf{r}^{\prime}_{j_{\ira}},\mathbf{r}^{\prime}_{j_{\ira}+1},\ldots,\mathbf{r}^{\prime}_{j_{\ira}+\bar{M}-1},\mathbf{r}^{\prime}_{j_{\ira}+\bar{M}})$
from the head of the worm to the bead $\alpha$ with the levy construction
where $\mathbf{r}^{\prime}_{j_{\ira}}=\mathbf{r}_{\ira}$ and
$\mathbf{r}^{\prime}_{j_{\ira}+\bar{M}}=\mathbf{r}_{\alpha}$.

The proposal probabilities by following such a scheme give

$\displaystyle\mathcal{P}_{\mathrm{Swap}}(\config\to\config^{\prime})$
$\displaystyle=\frac{\rho_{0}\left(\mathbf{r}_{\ira},\mathbf{r}_{\alpha};\bar{M}\tau\right)}{\Sigma_{1}}\frac{\prod_{k=0}^{\bar{M}-1}\rho_{0}\left(\mathbf{r}_{j_{\ira}+k}^{\prime},\mathbf{r}_{j_{\ira}+k+1}^{\prime};\tau\right)}{\rho_{0}\left(\mathbf{r}_{\ira},\mathbf{r}_{\alpha};\bar{M}\tau\right)}=\frac{\prod_{k=0}^{\bar{M}-1}\rho_{0}\left(\mathbf{r}_{j_{\ira}+k}^{\prime},\mathbf{r}_{j_{\ira}+k+1}^{\prime};\tau\right)}{\Sigma_{1}},$
(67) $\displaystyle\mathcal{P}_{\mathrm{Swap}}(\config^{\prime}\to\config)$
$\displaystyle=\frac{\rho_{0}\left(\mathbf{r}_{\ira},\mathbf{r}_{\alpha};\bar{M}\tau\right)}{\Sigma_{2}}\frac{\prod_{k=0}^{\bar{M}-1}\rho_{0}\left(\mathbf{r}_{j_{\ira}+k},\mathbf{r}_{j_{\ira}+k+1};\tau\right)}{\rho_{0}\left(\mathbf{r}_{\ira},\mathbf{r}_{\alpha};\bar{M}\tau\right)}=\frac{\prod_{k=0}^{\bar{M}-1}\rho_{0}\left(\mathbf{r}_{j_{\ira}+k},\mathbf{r}_{j_{\ira}+k+1};\tau\right)}{\Sigma_{2}},$
(68)

which results in the acceptance probabilities

$\displaystyle A_{\mathrm{Swap}}(\config\to\config^{\prime})$
$\displaystyle=\min\left(1,\frac{\Sigma}{\Sigma^{\prime}}e^{\tau\Delta
S(\config\to\config^{\prime})}\right),$ (69) $\displaystyle
A_{\mathrm{Swap}}(\config^{\prime}\to\config)$
$\displaystyle=\min\left(1,\frac{\Sigma^{\prime}}{\Sigma}e^{\tau\Delta
S(\config^{\prime}\to\config)}\right).$ (70)

Note, that we again perform importance sampling by choosing a worldline to
swap the worm with from the list $\chi$, thereby ensuring that the distance is
smaller than the order $\sqrt{\bar{M}\tau}$. This makes the update local,
which together with the fact that only a single permutation is made at a time
and modifications can be made many times before the configuration is closed,
ensures that exchanges or permutations are made if the system wants so.

Figure 10: Illustration of the swap update. The unchanged beads are portrayed
in black, the old beads and links are dashed, and the newly generated path or
Brownian bridge in solid grey.

The Markov chain updates, described above, combined with the canonical updates
(see section 1), form an ergodisic set of updates where each individual update
or complementary pair holds the detailed balance condition. This set allows
for conversion between the Z- and G-sector in a equilibriated and local
manner.

##### Worms algorithm in the Canonical ensemble

The set of updates in the Worm algorithm are in a natural way design to sample
the grand-canonical ensemble. Nevertheless, often one want to only take
measurement of a specific canonical ensemble $\\{N,V,T\\}$. One way of doing
this is choosing a optimal chemical potential $\mu$ such that the number of
particles fluctuates around the desired $N$ and only take measurements when
the system is in ensemble $\\{N,V,T\\}$171717This process is described in
section 4.. Nevertheless, this can be quite inefficient. Furthermore, this can
not be done if want the simulate non-interacting bosons as any non-zero
chemical potential will lead to a infinite of particles.

Another alternative is restricting the number of closed worldlines or particle
numbers to remain constant in the Z-sector. To enforce this on the Worm
algorithm, one often adds restriction on the set of updates. For example, if
one:

  * •

limits the set of the Worm algorithm to only include $\\{$Open, Close,
Advance, Recede, Swap$\\}$,

  * •

impede the time slices of head and tail to advance/recede passed each other,

  * •

confine $\bar{M}<J$,

one cannot change the particle number $N$ for configuration in the Z-sector.
Another option — which can lead to better performance — would be to design new
set of updates for the Worm algorithm. This is for example done in RotaPIMC.

###  4 Measure Observables

In the previous section, we explained how to generate new configurations in
the Markov chain. The only thing that remains is to extract physical
properties out of the configurations. As shown in section 1, this comes down
to determining the estimator $\mathcal{Q}(\config)$ of a thermodynamic
observable $\hat{O}$. This is done by writing out the definition (1)
explicitly using the small-time propagator you will be using in simulating the
system and, afterwards, comparing it to the form of equation 5. Often one has
also the option to extract it directly from the partition function by using
their thermodynamic definition, see for example section 1. If the estimator of
an observable is identified, it comes down to taking the average of this
estimator for all entries in the Markov chain as in equation 7.

We will cover the observables interesting for ultracold atoms in optical
lattices. However, note that WAPIMC is not limited by the observables covered
in this manuscript.

####  1 Energy

The estimator for the total energy of the system is often computed by using
the thermodynamic definition given by

$\displaystyle E=\expval{\ham}_{Z}=-\frac{1}{Z}\pdv{Z}{\beta}$

The subscript $Z$ emphasises the fact that the average is only taken over the
Z-sector. How this is explicitly done is described in section 4. The form of
the estimator is, of course, system dependent and hence relies on the content
of the density matrix. Let us compute the total energy estimator
$\mathcal{Q}_{E}$ for a non-interacting Bose gas in an external potential
described by the density matrix in (31)181818Note that we used the second-
order trotter formula to derive the density matrix and hence the result may
change if higher orders are used.. The derivation is analogous when adding
interaction to the system. The partition function can then be written as

$\displaystyle Z$
$\displaystyle=\int\dd{\config}\pi(\config)=\frac{1}{N!}\int\dd{\mathbf{R}_{0}}\ldots\dd{\mathbf{R}_{J-1}}\prod_{i=0}^{J-1}\rho_{0}\left(\mathbf{R}_{j},\mathbf{R}_{j+1};\tau\right)\mathrm{e}^{-S\left(\mathbf{R}_{j},\mathbf{R}_{j+1};\tau\right)}$
(71)
$\displaystyle=\frac{1}{N!}\int\dd{\mathbf{R}_{0}}\ldots\dd{\mathbf{R}_{J-1}}\frac{1}{\lambda_{\beta}^{dNJ}}\prod_{i=0}^{J-1}\exp\left(-\frac{\pi}{\lambda_{\tau}^{2}}\sum_{k=1}^{N}\left(\mathbf{r}_{k,j}-\mathbf{r}_{k,j+1}\right)^{2}-\frac{\tau}{2}\sum_{k=1}^{N}(V\left(\mathbf{r}_{k,j}\right)+V\left(\mathbf{r}_{k,j+1}\right))\right).$
(72)

where we remember that $\lambda_{\tau}=\sqrt{4\pi\lambda\tau}$ and $S$ is the
action. Substituting this partition function (71) in the expression for the
expectation value of the total energy (1) yields

$\displaystyle E=\frac{1}{Z}\int\dd{\config}\pdv{\pi(\config)}{\beta}.$ (73)

Further, one easily derives the identity

$\displaystyle\partial_{\beta}\left(\frac{1}{\lambda_{\beta}^{dNJ}}\right)=-\frac{dNJ}{2\beta}\left(\frac{1}{\lambda_{\beta}}\right)^{dNJ},$
(74)

combined with applying the product rule of derivatives to equation 73 we find
that the expectation value for the total energy of the system is given by

$\displaystyle E=$
$\displaystyle\frac{1}{Z}\int\dd{\config}\pi(\config)\times\left(\frac{dN}{2\tau}-\frac{1}{4\lambda\tau^{2}J}\sum_{j=0}^{J-1}\sum_{k=1}^{N}\left(\mathbf{r}_{k,j}-\mathbf{r}_{k,j+1}\right)^{2}+\frac{1}{2J}\sum_{j=0}^{J-1}\sum_{k=1}^{N}(V\left(\mathbf{r}_{k,j}\right)+V\left(\mathbf{r}_{k,j+1}\right))\right).$
(75)

Comparing this result to the form of equation 5, the total energy estimator
yields

$\displaystyle\mathcal{Q}_{E_{T}}=\frac{dN}{2\tau}-\frac{1}{4\lambda\tau^{2}J}\sum_{j=0}^{J-1}\sum_{k=1}^{N}\left(\mathbf{r}_{k,j}-\mathbf{r}_{k,j+1}\right)^{2}+\frac{1}{J}\sum_{j=0}^{J-1}\sum_{k=1}^{N}(V\left(\mathbf{r}_{k,j}\right)+V\left(\mathbf{r}_{k,j+1}\right)).$
(76)

The result is interpreted very easily. The first term corresponds to the well-
known _equipartition theorem_ for the kinetic energy of an ideal gas, e.g.,
$E_{\mathrm{kin}}=3Nk_{B}T/2$ in a three-dimensional system. The second term
is a quantum correction to kinetic energy. Indeed, in the classical limit
where only $J=1$ the correction term vanishes. The last term corresponds to
the contribution of potential energy. Notice that the kinetic energy becomes
small when the beads are far away from each other. It shows that worldlines in
imaginary time that combine to form a configuration do not represent real
physics, i.e., do not correspond to any real-time propagation, and are only is
a convenient representation of the quantum mechanical path integral structure.

Figure 11: Total energy of a boson in a harmonic trap
$V(\vb{r})=\frac{1}{2}\hbar\omega\vb{r}^{2}$ at different temperatures $T$ and
time slices $J$. The black dashed lines represent the exact value. The
coloured markers are the simulated values with the help of our implementation
of WAPIMC. Notice that the error bar becomes larger for increasing $J$.

The above estimator is known as the _thermodynamic estimator_ for the total
energy191919One could also compute the total energy estimator from the
explicit definition in equation 5 which will yield a different result for the
potential contribution. Nevertheless, both will converge to the true
expectation value in the limit $\tau\to 0$. A third estimator can also be
derived from the virial theorem. The combination of the three estimates can be
used to reduce the statistical error for the energy estimation Ceperley1995..
However, it is notorious to have a large error for a large $J$. This is due to
the the quantum correction term where in increasing $J$ the variance of the
sum diverges linearly Wolfhard1997. Indeed, in figure 11, we observe
empirically that the error for the total energy diverges when one uses time
slices for a boson in a harmonic potential. To circumvent the problem, we need
to propose a better kinetic energy estimator with a lower fluctuation. To this
goal, one often uses the _virial theorem_ ,
$\expval{\frac{\hat{p}^{2}}{2m}}=\frac{1}{2}\expval{\hat{x}V^{\prime}(\hat{x})}$
where we have used the abbreviation $V^{\prime}(x)=dV/dx$. Using the
expression, one can derive a “variance reduced” or _virial estimator_ , given
by

$\displaystyle\mathcal{Q}_{E_{V}}=\frac{dN}{2\beta}$
$\displaystyle+\frac{1}{4\lambda\tau
J}\sum_{k=1}^{N}\left(\mathbf{r}_{k,J-1}-\mathbf{r}_{k,J}\right)\left(\mathbf{r}_{k,J}-\mathbf{r}_{k,0}\right)$
(77)
$\displaystyle+\frac{1}{2\beta}\sum_{j=1}^{J-1}\sum_{k=1}^{N}\left(\mathbf{r}_{k,j}-\mathbf{r}_{k,0}\right)\frac{\partial}{\partial\mathbf{r}_{k,j}}\left(S\left(\mathbf{r}_{k,j+1},\mathbf{r}_{k,j};\tau\right)+S\left(\mathbf{r}_{k,j},\mathbf{r}_{k,j-1};\tau\right)\right)$
(78) $\displaystyle+\sum_{j=0}^{J-1}\sum_{k=1}^{N}\frac{\partial
S\left(\mathbf{r}_{k,j+1},\mathbf{r}_{k,j};\tau\right)}{\partial\tau},$ (79)

with $S$ the action or potential contribution to the density matrix and the
bead notation is always understood the $\beta$-periodic sense. This expression
has been taken from Spada2022. Although, different versions are achievable,
e.g., in RotaPIMC. For a derivation and comparison between different energy
estimators, Ceperley1995 is recommended.

We conclude with two final remarks. First, the second term in the virial
estimator (77) represents the particle statistics. Indeed, if no permutations
between the worldlines are present, we have that
$\mathbf{r}_{k,J}=\mathbf{r}_{k,0}$ and the term vanishes. Next, the
difference $\mathbf{r}_{k,J}-\mathbf{r}_{k,0}$ in the second must be evaluated
continuously along the worldline202020Periodic boundary conditions may only
apply between two succeeding beads!. It follows, that the absolute value of
this difference can be larger than the size of the simulation box.

####  2 Spatial structure

##### Density

The spatial density distribution $n(\mathbf{r})$ of a system often reveals
interesting features of the system. It turns out the density of the system is
one of the variables that is easily obtained as it is diagonal in the position
representation. This allows us to easily compute the definition of the
estimator (5) to be
$\mel{\vb{R}_{0}}{\hat{n}(\vb{r})}{\vb{R}_{0}}=\sum^{N}_{i=1}\delta(\vb{r}-\vb{r}_{i,0})$.
However, we can improve the statistics error by noticing that the translation
invariance in imaginary time makes that a time slice is not particularly
special. As a result, an improved version is found for the local spatial
density estimator at a position $\mathbf{r}_{\alpha}$, given by

$\displaystyle\mathcal{Q}_{n}(\mathbf{r}_{\alpha})=\frac{1}{J}\sum_{i=1}^{N}\sum_{j=0}^{J-1}\frac{\delta\left(\vb{r}_{\alpha},\mathbf{r}_{i,j}\right)}{V_{\alpha}}.$
(80)

Here, the spatial volume is divided in a discrete grid of hyper-cubes in $d$
dimensions where each bin $\alpha$ has a center coordinate
$\mathbf{r}_{\alpha}$ with volume $V_{\alpha}$ such that

$\displaystyle\delta\left(\mathbf{r}_{j,k},\mathbf{r}_{\alpha}\right)=\begin{cases}1,&\qq{if}\mathbf{r}_{i,j}\in\alpha,\\\
0,&\text{ otherwise }.\end{cases}$ (81)

##### Radial Distribution

Using the same approach as for the density, other spatial observables can
easily be extracted from the configuration, e.g., the radial distribution
function $g(r)$, a pair correlation function that explains how the atoms in a
system are radially distributed around each other. This is a valuable way of
representing the average structure of disordered systems as it can be easily
measured experimentally. More concretely, it is proportional to the
probability of finding a particle at a distance of $r$ away from a given
reference particle, which yields

$\displaystyle\mathcal{Q}_{g}({r}_{\alpha})=\frac{V}{N^{2}J}\sum_{j=0}^{J-1}\sum_{i\neq
k}\frac{\delta\left({r}_{\alpha}-\left({r}_{i,j}-{r}_{k,j}\right)\right)}{L_{\alpha}},$
(82)

evaluating the delta-function by $1$ if ${r}_{\alpha}$ lies in the interval
$\alpha=[r_{\alpha}-\Delta r,r_{\alpha}+\Delta r]$ with length
$L_{\alpha}=2\Delta r$. In general, one could look many other pair correlation
function such as the Static structure factor, etc RotaPIMC.

##### Compressibility

The spatial density $n$ in combination with working in the grand canonical
ensemble gives us the opportunity to compute the compressibility $\kappa$ of
the system. It can be defined as the change in spatial density $n$ as the
chemical potential $\mu$ is varied at a constant temperature, i.e.,

$\displaystyle\kappa=\left.\frac{\partial n}{\partial\mu}\right|_{T}.$ (83)

Substituting the grand canonical spatial density operator
$\hat{n}=\operatorname{Tr}\left[\hat{N}e^{-\beta(\hat{H}-\mu\hat{N})}\right]/VZ$
yields that the compressibility $\kappa$ is dependent on the variance of the
number of particles, i.e.,Yao2020

$\displaystyle\kappa=V\beta\frac{\expval{N^{2}}_{\mathrm{Z}}-\expval{N}_{\mathrm{Z}}^{2}}{\expval{N}_{\mathrm{Z}}}.$
(84)

Here the subscript $Z$ indicates the average in equation 7 only includes a
configuration in the Z-sector. The estimators $\mathcal{Q}_{N}$ and
$\mathcal{Q}_{N^{2}}$ in (7) are just the number of closed worldlines and the
number of closed worldlines squared in a configuration $\config$.

####  3 Superfluid density

As PIMC solves the path integral formulation of quantum mechanics also pure
quantum mechanical properties can be predicted. An example would be
superfluidity. Nevertheless, there is at the moment no complete satisfying
coherent description theoretical description of superfluidity212121The second
chapter of Boning2007 contains a nice overview of the challenges.. It is
currently well understood in a phenomenological context, although analytical
conclusions obtained from microscopic models remain sporadic Boning2007. For
our purposes, it is enough to understand superfluidity as the decrease of the
system’s moment of inertia induced by quantum effects, i.e. non-classical
rotational inertia. In what follows, we will derive an estimator both for the
global superfluid fraction and local superfluid density with the help of
Ceperley1995; Boning2007; Dornheim2018.

##### Global superfluid fraction

Consider the system is placed in a bucket and is slowly rotated with an
angular frequency $\omega$ around the z-axis. The bucket does not have to be
cylindrical symmetric. A classical fluid will rotate together with the walls
where the work done is $\frac{1}{2}I_{\mathrm{c}}\omega^{2}$ with the
$I_{\mathrm{c}}$ the classical moment of inertia. The cycle duration is
thereby proportional to the square root of the moment of inertia and,
therefore, to the total density n, i.e.

$\displaystyle\frac{2\pi}{\omega}\propto\sqrt{I_{\mathrm{c}}}\propto\sqrt{n}.$
(85)

Assume that at low temperature, only a portion of the density, dubbed
$n_{\mathrm{n}}$, will participate with the rotation, such that

$\displaystyle\sqrt{I(T)}\propto\sqrt{n_{\mathrm{n}}}=\sqrt{n-n_{\mathrm{sf}}},$
(86)

with $n_{\mathrm{sf}}$ the density of the fluid which does not respond to the
rotation. What follows is the relation

$\displaystyle
f\equiv\frac{n_{\mathrm{sf}}}{n}=1-\frac{I(T)}{I_{\mathrm{c}}},$ (87)

where we defined the superfluid fraction $f$ as the ratio of the superfluid
and total density. What remains is to find an expression for the momenta of
inertia. Following Sindzingre1989 and Ceperley1987 one can use the definition

$\displaystyle
I=\left.\dv{F}{\omega^{2}}\right|_{\omega=0}=\left.\dv{\left\langle
L_{\mathrm{z}}\right\rangle}{\omega}\right|_{\omega=0}$

where the momenta of inertia is the work done for an infinitesimally small
rotation rate. Rewriting the expectation value of the angular momentum
operator $L_{\mathrm{z}}$ in the path integral formulation and decomposition
in $J$ time slices yields the following expression Ceperley1987:

$\displaystyle\mathcal{Q}_{\mathrm{f}}=\frac{2m}{\lambda\beta}\frac{\left\langle
A_{z}^{2}\right\rangle}{I_{\mathrm{c}}},$ (88)

with

$\displaystyle\mathbf{A}=\frac{1}{2}\sum_{i=1}^{N}\sum_{j=0}^{J-1}\left(\mathbf{r}_{i,j}\times\mathbf{r}_{i,j+1}\right)\qq{and}I_{c}=m\left\langle\sum_{i=1}^{N}\sum_{j=0}^{J-1}\mathbf{r}_{i,j}^{\perp}\cdot\mathbf{r}_{i,j+1}^{\perp}\right\rangle.$
(89)

The vector $\vb{A}$ can be seen as the area of the path or ring-polymers.
Hence, equation (88) connects the superfluid fraction to the worldlines mean-
squared area normalized by the cross-sectional area of the whole system. This
links nicely to the particle exchange picture as the permutated polymers can
cover a larger area. Taking a system where the bosons can only move in two
spatial dimensions, we can write the superfluid fraction as

$\displaystyle\mathcal{Q}_{\mathrm{f}}=\frac{1}{\lambda\beta}\frac{\sum_{i=1}^{N}\sum_{j=0}^{J-1}r^{x}_{i,j}r^{y}_{i,j+1}-r^{y}_{i,j}r^{x}_{i,j+1}}{\sum_{i=1}^{N}\sum_{j=0}^{J-1}r^{x}_{i,j}r^{x}_{i,j+1}+r^{y}_{i,j}r^{y}_{i,j+1}},$
(90)

where the rotation axis has been taken as perpendicular to the two-dimensional
plane. Note that the superfluid fraction is indeed dimensionless.

##### Local superfluid density

It is useful to have an estimator that measures the global properties of the
system. However, acquiring spatially resolved information about the system of
interest is often also necessary. For the superfluid density, this is provided
by Kwon2006. Although previous efforts were made by Draeger2003, Kwon2006
showed their implementation leads to biased results as the integration does
not give the correct moment of inertia $I$.

The estimator works on the basis of the Landau two fluid model, which divides
total density into normal $n_{\mathrm{n}}$ and superfluid contributions of
$n_{\mathrm{sf}}$, i.e., $n=n_{\mathrm{n}}+n_{\mathrm{sf}}$. Therefore, we can
write the moment of inertia as

$\displaystyle
I(T)=m\int\dd{\mathbf{r}}\left(n(\mathbf{r})-n_{\mathrm{sf}}(\mathbf{r})\right)\mathbf{r}_{\perp}^{2}=I_{\mathrm{cl}}-m\int\dd{\mathbf{r}}n_{\mathrm{sf}}(\mathbf{r})\mathbf{r}_{\perp}^{2},$
(91)

where we used that the classical moment of inertia
$I_{\mathrm{cl}}=m\int\dd{\mathbf{r}}n(\mathbf{r})\mathbf{r}_{\perp}^{2}$.
$\mathbf{r}_{\perp}$ is the vector where we integrate over with the components
of the rotation axis (the z-axis) set to zero. Using the relation (87) and the
estimator (88) that we deduced in the previous subsection, we can write

$\displaystyle
m\int\dd{\mathbf{r}}n_{\mathrm{sf}}(\mathbf{r})\mathbf{r}_{\perp}^{2}=\frac{2m}{\lambda\beta}\left\langle
A_{z}\right\rangle^{2}.$ (92)

with $A_{z}$ the z-component of the area vector defined in (89). Above
relation yields

$\displaystyle\mathcal{Q}_{n_{\mathrm{sf}}}(\mathbf{r})=\frac{2}{\lambda\beta}\frac{A_{z}A_{z,\mathrm{loc}}(\mathbf{r})}{\vb{r}_{\perp}^{2}},\qq{with}\mathbf{A}_{\mathrm{loc}}(\mathbf{r})=\frac{1}{2}\sum_{i=1}^{N}\sum_{j=0}^{J-1}\left(\mathbf{r}_{i,j}\times\mathbf{r}_{i,j+1}\right)\delta\left(\mathbf{r}-\mathbf{r}_{i,j}\right).$
(93)

The ladder constitutes the local contribution of all trajectories, as it is
the sum of all area segments travelling through a differential volume enforced
by the delta function. The delta function is dealt with in the same manner as
the density, i.e., by dividing space into discrete hyper-cubes and using
equation (81).

####  4 Measuring scheme

Figure 12: Illustration of two different measuring schemes red(a) and red(b).
Each square in the two series describes an element of the Markov chain which
contains diagonal (green) and off-diagonal (red) configurations. The black
arrow corresponds to attempting a measurement. The concept of the figure is
based on dornheim2014.

As explained in section 1, the canonical expectation value of an observable
$\oper$ is given by

$\displaystyle\expval*{\oper}\approx\frac{\mathcal{Q}\left(\config_{1}\right)+\ldots+\mathcal{Q}\left(\config_{n_{\mathrm{MC}}}\right)}{n_{\mathrm{MC}}},$
(94)

where the set
$\\{\config_{k}\\}=\\{\config_{1},\ldots,\config_{N_{\mathrm{MC}}}\\}$
represents the the Markov chain of canonical configurations. However, WAPIMC
generates configuration in an extended configuration space. Hence, we need a
measuring scheme to extract the canonical date from the Markov chain
$\\{\config_{k}\\}$.

In figure 12, two different measuring schemes represented by red(a) and red(b)
are shown. In general both schemes are equally valid and which is preferred
depends on the system. Measure scheme red(a) goes as follows:

  * •

Keep track of the amount of elements in the Markov chain.

  * •

Check if the walker is in the z-sector (green square) every $k$ updates (black
arrow). If affirmative then perform the canonical meaurements and add them to
the corresponding canonical set. If the walker is in G-sector (red square)
then do nothing

  * •

Reiterate until enough measurements are taken for the desired canonical set of
measurements.

The second option would be:

  * •

Count every new diagonal configuration in the Markov chain with a specific
number of particle $N$ stored in $\operatorname{ctr}(N)$.

  * •

If one of the counters $\operatorname{ctr}(N)$ is divisible by $k$, then
perform the canonical measurements and add them to the corresponding canonical
set. Afterwards, reset the counter $\operatorname{ctr}(N)$.

  * •

Reiterate until enough measurements are taken for the desired canonical set of
measurements.

Both meaures schemes are equally valid. Which of the two extraction strategies
is more efficient depends on the specific use case. The first strategy is
applied in this work.

####  5 Systematic and Statistical errors in PIMC

As in every simulation, we need error bars in our simulation results to assess
if the obtained results are decent. The following summary of the error
analyses of a PIMC simulation is based on Vinay2010; Young2014; Janke2004. The
thesis’s dornheim2014; Yao2020 also provide a section on the topic.

The PIMC simulation results include two types of errors that are common. The
first is the systematic error caused by the inaccuracy of the density matrix,
an error that goes as an increasing function of $\tau$. This finite value has
an effect on any result that uses the small-time propagator. The error is
often handled by observing the convergence of the system for $\tau\to 0$ where
the result is extrapolated.

The other type of error is a statistical error by virtue of the limited number
of configurations generated during the simulation, contributing to the
expectation value of the observable in equation 7. As stated earlier in
section 2, if there are a sufficient enough number of uncorrelated
configurations, the central limit theorem states that the obtained average
$\bar{O}$ of an observable $\hat{O}$ is normally distributed around the exact
expectation value with a finite variance given by

$\displaystyle\sigma_{\bar{O}}^{2}=\expval{\bar{O}^{2}}-\expval{\bar{O}}^{2}=\frac{\sigma_{O_{i}}^{2}}{n},$
(95)

where $n$ is the number of Markov chain elements and $\sigma_{{O}_{k}}$ the
variance of the individual measurement ${O}_{i}$ with $i\in{1,\ldots,n}$. It
implies the notorious Monte Carlo property that the error vanishes as
$1/\sqrt{n}$.

##### Autocorrelation

Nevertheless, the result is different if we do not neglect the correlation
between the individual configurations. This is easily observed if we
substitute the expression for the mean $\bar{O}=\sum^{n}_{i=1}O_{i}/n$ into
equation 95 which yields

$\displaystyle{\sigma_{\bar{O}}}^{2}=\frac{1}{n^{2}}\sum_{i=1}^{n}\left(\expval{O_{i}^{2}}-\expval{O_{i}}^{2}\right)+\frac{1}{n^{2}}\sum_{i\neq
k}^{n}\left(\expval{O_{i}O_{k}}-\expval{O_{i}}\expval{O_{j}}\right).$ (96)

From the first term, we identify the variance $\sigma_{O_{i}}^{2}$ of an
individual Markov chain element. The result can be simplified further by using
the summation of the summation in the second term to write $\sum_{i\neq
k}^{n}=2\sum_{i=1}^{n}\sum_{k=i+1}^{n}$. Furthermore, using the time
translation invariance, we obtain

$\displaystyle{\sigma_{\bar{O}}}^{2}=\frac{1}{n}\left[\sigma_{O_{i}}^{2}+2\sum_{k=1}^{n}\left(\expval{O_{1}O_{1+k}}-\expval{O_{1}}\expval{O_{1+k}}\right)\left(1-\frac{k}{n}\right)\right].$
(97)

One can factor out $\sigma_{O_{i}}^{2}$ to write the variance of the mean in
its usual form, i.e.,

$\displaystyle{\sigma_{\bar{O}}}^{2}=2\frac{\sigma_{O_{i}}^{2}}{n}\tau_{\mathrm{auto}}$
(98)

, where the integrated autocorrelation time $\tau_{\mathrm{auto}}$ represents
a measure for the length between two independent configurations in the Markov
chain, given by

$\displaystyle\tau_{\mathrm{auto}}=\frac{1}{2}+\sum_{k=1}^{n}O(k)\left(1-\frac{k}{n}\right)\qq{with}O(k)=\frac{\expval{O_{i}O_{i+k}}-\expval{O_{i}}\expval{O_{i}}}{\expval{O_{i}^{2}}-\expval{O_{i}}\expval{O_{i}}}.$
(99)

It is clear that correlations between configurations increases the variance.