Imaginary-time nonuniform mesh method for solving the multidimensional Schr\"odinger equation: Fermionization and melting of quantum Lennard-Jones crystals
Alberto Hernando
An imaginary-time nonuniform mesh method is presented and used to find the first 50 eigenstates and energies of up to five strongly interacting spinless quantum Lennard-Jones particles trapped in a one-dimensional harmonic potential. We show that the use of tailored grids reduces drastically the computational effort needed to diagonalize the Hamiltonian and results in a favorable scaling with dimensionality. Solutions to both bosonic and fermionic counterparts of this strongly interacting system are obtained, the bosonic case clustering as a Tonks-Girardeau crystal exhibiting the phenomenon of fermionization. The numerically exact excited states are used to describe the melting of this crystal at finite temperature.
tise_method_v11.tex
\documentclass[prl,twocolumn]{revtex4-1}%
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage[latin9]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[english]{babel}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{color}%
\setcounter{MaxMatrixCols}{30}
%TCIDATA{OutputFilter=latex2.dll}
%TCIDATA{Version=5.50.0.2960}
%TCIDATA{LastRevised=Tuesday, June 04, 2013 20:01:27}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="SaveForMode" CONTENT="2">}
%TCIDATA{BibliographyScheme=BibTeX}
%TCIDATA{Language=American English}
%BeginMSIPreambleData
\providecommand{\U}[1]{\protect\rule{.1in}{.1in}}
%EndMSIPreambleData
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\ben}{\begin{eqnarray}}
\newcommand{\een}{\end{eqnarray}}
\newcommand{\q}{\ensuremath{{\mathbf q}}}
%BeginMSIPreambleData
\ifx\pdfoutput\relax\let\pdfoutput=\undefined\fi
\newcount\msipdfoutput
\ifx\pdfoutput\undefined\else
\ifcase\pdfoutput\else
\msipdfoutput=1
\ifx\paperwidth\undefined\else
\ifdim\paperheight=0pt\relax\else\pdfpageheight\paperheight\fi
\ifdim\paperwidth=0pt\relax\else\pdfpagewidth\paperwidth\fi
\fi\fi\fi
%EndMSIPreambleData
\begin{document}
\title{Imaginary-time nonuniform mesh method for solving the multidimensional
Schr\"odinger equation: Fermionization and melting of quantum Lennard-Jones crystals}
\author{Alberto Hernando}
\email{alberto.hernandodecastro@epfl.ch}
\author{Ji\v{r}\' i Van\'i\v{c}ek}
\email{jiri.vanicek@epfl.ch}
\affiliation{Laboratory of Theoretical Physical Chemistry, Institut des Sciences et
Ing\'enierie Chimiques, \'Ecole Polytechnique F\'ed\'erale de Lausanne,
CH-1015 Lausanne, Switzerland}

\begin{abstract}
An imaginary-time nonuniform mesh method is presented and used to find the
first 50 eigenstates and energies of up to five strongly interacting spinless
quantum Lennard-Jones particles trapped in a one-dimensional harmonic
potential. We show that the use of tailored grids reduces drastically the
computational effort needed to diagonalize the Hamiltonian and results in a
favorable scaling with dimensionality. Solutions to both bosonic and fermionic
counterparts of this strongly interacting system are obtained, the bosonic
case clustering as a Tonks-Girardeau crystal exhibiting the phenomenon of
fermionization. The numerically exact excited states are used to describe the
melting of this crystal at finite temperature.

\end{abstract}
\date{\today}
\maketitle


The multidimensional Schr\"odinger equation (MDSE) is undoubtedly one of the
cornerstones of modern physics and much attention has been paid to developing
efficient numerical methods for finding its solutions
\cite{Harris:1965,Kosloff:1988,Shimshovitz:2012,Pollet:2012,book_Parr_Yang,
Dalfovo:1995,Savenko:2013,Rogel-Salazar:2013,Billing:2000,Runge_Gross:1984,Abid:2003,book_Wyatt,
Alvarez:2011,Tkatchenko:2012,Delgado:1996,Corney:2004,Torres-Vega:1991,McMahon:2012}%
. A very rich testing ground for such methods has been provided by the
observation of new quantum phases at ultracold temperatures in finite and
homogeneous systems
\cite{McMahon:2012,Kapitza:1938,Leggett:1972,Grebenev:2000,Anglin_Ketterle:2002,Balibar:2010}%
, and also by the development of optical lattices where ultracold atoms are
trapped \cite{Bloch:2005}. Due to their fascinating structural and dynamical
properties, special attention has been recently devoted to one-dimensional
traps
\cite{Ronzheimer:2013,Panfil:2013,Vignolo:2013,Gring:2012,Kinoshita:2006,Paredes:2004}%
. Indeed, in the strongly interacting (Tonks-Girardeau) regime of bosonic
particles trapped in one-dimensional geometries, the repulsive nature of the
atomic interaction at short distances gives rise to the phenomenon known as
\textit{fermionization}, the mechanism of which is actively studied both
theoretically and experimentally. 

Rigorous description and explanation of the new physics found in these
well-controlled experiments require accurate theoretical methods and
constitute a formidable challenge \cite{Bloch:2008}, the main technical
difficulty being the scaling of numerical algorithms with the number of
dimensions $D$. Indeed, standard algorithms for solving differential
equations, such as the Finite Difference method, scale exponentially with
dimensions \cite{book_Morton_Mayers}, making numerical solutions of
many-dimensional problems impracticable, if not impossible. Improved methods
addressing this difficulty in the case of stationary states include the
Discrete Variable Representation (DVR) \cite{Harris:1965}, collocation method
\cite{Kosloff:1988}, phase-space method based on von Neumann periodic lattice
\cite{Shimshovitz:2012}, variational or diffusion quantum Monte Carlo (MC)
methods \cite{McMahon:2012,Pollet:2012}, Density Functional Theory (DFT)
\cite{McMahon:2012,book_Parr_Yang}, mean-field or pseudopotential interaction
models \cite{Dalfovo:1995,Savenko:2013,Rogel-Salazar:2013}, and many others.
Some of these methods find only the ground state of the time-independent MDSE,
using different efficient techniques such as the imaginary time (IT)
propagation \cite{Popov:2005} or the Variational Principle \cite{book_Sakurai}%
. Methods for real-time quantum dynamics include the Time-Dependent DVR
\cite{Billing:2000}, DFT \cite{Runge_Gross:1984}, mean-field approaches
\cite{Abid:2003}, trajectory-based methods such as Bohmian dynamics
\cite{book_Wyatt}, or time-dependent density matrix renormalization group
(t-DMRG) method \cite{Alvarez:2011}, which has proven to be very efficient in
one-dimensional geometries. Despite many accomplishments in special cases,
finding excited states and describing the real-time dynamics governed by a
general high-dimensional Hamiltonian in the strongly interacting regime
remains a difficult computational challenge.

In this paper \textsl{we propose a novel general method, scaling favorably
with dimensions, which is able to solve the time-independent MDSE numerically
exactly and simultaneously finds both its ground and excited states}.
Obviously, the proposed IT nonuniform mesh method (ITNUMM) is not intended to
replace other well established approaches; instead we expect it to have a
domain of applicability where other methods present more technical
difficulties, such as in finding excited states of many-dimensional systems
and where efficiency is more important than high accuracy. To show that ITNUMM
achieves these goals, we apply it to find the wavefunctions of the first 50
states of an ensemble of up to five distinguishable Lennard-Jones (LJ)
spinless particles trapped in a one-dimensional harmonic potential in the
Tonks-Girardeau regime. Once these states are obtained, we find, via
symmetrization and anti-symmetrization, the solutions for the Bose-Einstein
and Fermi-Dirac statistics, respectively, and observe fermionization in the
bosonic case. We also show that the computed excited states can be used in a
thermal average to describe the melting of the LJ\ clusters at finite
temperature. \textsl{As we use no other approximation than the numerical
discretization of space and time, the obtained results are numerically exact}.

The derivation of our method starts by rewriting the time-dependent MDSE
\cite{book_Sakurai}
\begin{equation}
i\hbar\frac{d}{dt}|\psi(t)\rangle=\mathcal{H}|\psi(t)\rangle,
\end{equation}
where $|\psi(t)\rangle$ is the quantum state at time $t$ of the $D$%
-dimensional system described by Hamiltonian $\mathcal{H}$, in terms of the
quantum propagator $K({\mathbf{q}},{\mathbf{q}}^{\prime};t-t^{\prime
}):=\langle{\mathbf{q}}|e^{-i(t-t^{\prime})\mathcal{H}/\hbar}|{\mathbf{q}%
}^{\prime}\rangle$ in the position basis $|{\mathbf{q}}\rangle$:%
\begin{equation}
\psi({\mathbf{q}},t)=\int d{\mathbf{q}}^{\prime}K({\mathbf{q}},{\mathbf{q}%
}^{\prime};t-t^{\prime})\psi({\mathbf{q}}^{\prime},t^{\prime}).
\end{equation}
Hamiltonian $\mathcal{H}:=\mathcal{H}_{0}+\mathcal{H}_{1}$ is now split into
two components: $\mathcal{H}_{0}$ is any Hamiltonian that includes the kinetic
energy operator $\mathcal{T}$ and whose matrix elements in the ${\mathbf{q}}%
$-representation are known, while $\mathcal{H}_{1}\equiv\mathcal{H}%
_{1}({\mathbf{q}})$ is any many-body potential depending only on ${\mathbf{q}%
}$. For very short time intervals $t-t^{\prime}=\Delta t$, the time evolution
operator can be split to first order as $e^{-i\Delta t\mathcal{H}/\hbar
}=e^{-i\Delta t\mathcal{H}_{0}/\hbar}e^{-i\Delta t\mathcal{H}_{1}/\hbar
}+O(\Delta t^{2})$ and one can write
\begin{align}
\psi({\mathbf{q}},t^{\prime}+\Delta t)  &  = \int d{\mathbf{q}}^{\prime}%
K_{0}({\mathbf{q}},{\mathbf{q}}^{\prime};\Delta t)e^{-i\Delta t\mathcal{H}%
_{1}({\mathbf{q}}^{\prime})/\hbar}\psi({\mathbf{q}}^{\prime},t^{\prime
})\nonumber\\
&  \ +O(\Delta t^{2}), \label{eq:1}%
\end{align}
where $K_{0}({\mathbf{q}},{\mathbf{q}}^{\prime};\Delta t):=\langle{\mathbf{q}%
}|e^{-i\Delta t\mathcal{H}_{0}/\hbar}|{\mathbf{q}}^{\prime}\rangle$ is the
propagator of $\mathcal{H}_{0}$, which is assumed to be known explicitly.

The $|{\mathbf{q}}\rangle$ basis is discretized as
\begin{equation}
\int d{\mathbf{q}}|{\mathbf{q}}\rangle\langle{\mathbf{q}}|=\lim_{N\rightarrow
\infty}\sum_{j=1}^{N}w({\mathbf{q}}_{j})|{\mathbf{q}}_{j}\rangle
\langle{\mathbf{q}}_{j}|
\end{equation}
where $w({\mathbf{q}}_{j})$ is a weight function depending on a particular
realization of the $N$ states $|{\mathbf{q}}_{j}\rangle$. Indeed, $w$ is
defined as $w({\mathbf{q}}):=[Np({\mathbf{q}})]^{-1}$, where $p({\mathbf{q}})$
is the density distribution of the ${\mathbf{q}}_{j}$. With this
discretization, Eq. (\ref{eq:1}) becomes
\begin{align}
&  \!\!\!\!\psi({\mathbf{q}}_{j},t^{\prime}+\Delta t)=\label{eq:2}\\
&  \lim_{N\rightarrow\infty}\sum_{k=1}^{N}w({\mathbf{q}}_{k})K_{0}%
({\mathbf{q}}_{j},{\mathbf{q}}_{k};\Delta t)e^{-i\Delta t\mathcal{H}%
_{1}({\mathbf{q}}_{k})/\hbar}\psi({\mathbf{q}}_{k},t^{\prime}).\nonumber
\end{align}
Since our main interest is finding the stationary states of $\mathcal{H}$, in
the following we will assume that (i) $\psi({\mathbf{q}},t)=e^{-itE_{n}/\hbar
}\varphi_{n}({\mathbf{q}})$ where $\varphi_{n}({\mathbf{q}})$ and $E_{n}$ are
the \mbox{$n$th} eigenstate and eigenenergy of the Hamiltonian $\mathcal{H}$,
and that (ii) the evolution is performed in IT ($t\rightarrow-i\tau$).
Although the density $p({\mathbf{q}})$ is arbitrary, below we show that Eq.
(\ref{eq:2}) simplifies in the IT scheme if this density corresponds to the
classical Boltzmann distribution of $\mathcal{H}_{1}$, namely if
\begin{equation}
p({\mathbf{q}})=Z_{\mathcal{H}_{1}}^{-1}e^{-\Delta\tau\mathcal{H}%
_{1}({\mathbf{q}})/\hbar},
\end{equation}
where $Z_{\mathcal{H}_{1}}=\mathrm{Tr}e^{-\Delta\tau\mathcal{H}_{1}/\hbar}$ is
a normalization constant (called configuration integral) and $\Delta\tau
/\hbar$ plays the role of the inverse temperature $\beta$. Under these
conditions, Eq. (\ref{eq:2}) reads
\begin{align}
&  \!\!\!\!\!\!\!\!e^{-\Delta\tau E_{n}/\hbar}\varphi_{n}({\mathbf{q}}%
_{j})=\label{eq:3}\\
&  \lim_{N\rightarrow\infty}\frac{Z_{\mathcal{H}_{1}}}{N}\sum_{k=1}^{N}%
K_{0}({\mathbf{q}}_{j},{\mathbf{q}}_{k};-i\Delta\tau)\varphi_{n}({\mathbf{q}%
}_{k}).\nonumber
\end{align}
By defining vector $\Phi_{n}:=\{\varphi_{n}({\mathbf{q}}_{j})\}_{j=1}^{N}$,
whose $j$th component is the wavefunction evaluated at position ${\mathbf{q}%
}_{j}$, and matrix $\hat{K}_{jk}:=K_{0}({\mathbf{q}}_{j},{\mathbf{q}}%
_{k};-i\Delta\tau)Z_{\mathcal{H}_{1}}/N$ whose elements are proportional to
the propagator $K_{0}$ from ${\mathbf{q}}_{j}$ to ${\mathbf{q}}_{k}$, one can
rewrite Eq. (\ref{eq:3}) as a matrix eigenvalue equation
\begin{equation}
e^{-\Delta\tau E_{n}/\hbar}\Phi_{n}=\hat{K}\cdot\Phi_{n}.
\end{equation}
This equation, central to the ITNUMM, exhibits the main advantage of our
method---\textsl{the problem of finding the spectrum and eigenfunctions of the
original Hamiltonian $\mathcal{H}$ is reduced to sampling the classical
Boltzmann distribution and diagonalizing $\hat{K}$ evaluated at those points}.
Instead of the Hamiltonian, we diagonalize the imaginary-time propagator,
i.e., a matrix with analytically known and real-valued elements.
\textsl{Evaluation of, e.g., derivatives or Fourier transforms is not needed}.
Indeed, the implementation of the algorithm is rather simple since it only
requires standard methods for sampling from arbitrary probability
distributions and diagonalizing sparse real-valued matrices. The computational
effort is also reduced by constructing a nonuniform grid in which more grid
points are placed in areas where the wavefunctions exhibit more detailed
features. In the special case of $\mathcal{H}_{0}\equiv\mathcal{T}$,
$\mathcal{H}_{1}({\mathbf{q}})$ equals the classical potential energy, $K_{0}$
is a free-particle propagator in $D$ dimensions \cite{book_Sakurai}, and
matrix elements $\hat{K}_{jk}$ assume the Gaussian form
\begin{equation}
\hat{K}_{jk}=\frac{Z_{\mathcal{H}_{1}}}{N}\left(  \frac{m}{2\pi\hbar\Delta
\tau}\right)  ^{D/2}\exp\left[  -\frac{m}{2\hbar\Delta\tau}({\mathbf{q}}%
_{j}-{\mathbf{q}}_{k})^{2}\right]  ,
\end{equation}
where $m$ is the mass, for simplicity assumed to be the same for all degrees
of freedom. In correlated systems, where sampling the Boltzmann distribution
is difficult or unfeasible---as in the case of Coulomb interaction, we propose
the splitting $\mathcal{H}_{1}({\mathbf{q}})=V_{1}({\mathbf{q}})+V_{2}%
({\mathbf{q}})$, where $V_{1}({\mathbf{q}})$ is a sum of well-behaved one-body
potentials and $V_{2}({\mathbf{q}})$ is the remainder including all
correlations. Here the sampling is performed with the weight $p({\mathbf{q}%
})=Z_{V_{1}}^{-1}e^{-\Delta\tau V_{1}({\mathbf{q}})/\hbar}$ and normalization
$Z_{V_{1}}=\mathrm{Tr}e^{-\Delta\tau V_{1}/\hbar}$; the matrix to be
diagonalized becomes
\begin{align}
\hat{K}_{jk}  &  =\frac{Z_{V_{1}}}{N}\left(  \frac{m}{2\pi\hbar\Delta\tau
}\right)  ^{D/2}\nonumber\\
&  \times\exp\left[  -\frac{m}{2\hbar\Delta\tau}({\mathbf{q}}_{j}-{\mathbf{q}%
}_{k})^{2}-\frac{\Delta\tau}{\hbar}V_{2}({\mathbf{q}}_{k})\right]  .
\label{eq:4}%
\end{align}


We have found this method to be very efficient in one-dimensional problems
with several very different potentials. Although an arbitrary sampling
procedure can be used, we have employed a \textit{quadrature scheme}: instead
of random sampling of $p(q)$ by a MC procedure, the $q_{j}$ points are chosen
with a deterministic algorithm. The motivation for this approach is reducing
to a minimum the number of vector-elements needed for a given accuracy, and
thus reducing the computational cost of the diagonalization of $\hat{K}$.
Specifically, we first consider a new variable $u$, uniformly distributed in
the interval $[0,1]$, and define an equidistant grid $u_{j}=(j-1/2)/N$. The
Jacobian of the transformation from $q$ to $u$ is given by $p(q)$ since
$p(q)dq=du$, hence
\begin{equation}
u(q)=\int_{-\infty}^{q}dq^{\prime}p(q^{\prime})=P(q),
\end{equation}
where $P(q)$ is the cumulative distribution function. Next, the $q$-grid is
obtained by inverting this equation for all values of $u_{j}$, and once the
$q$-grid is ready, the evaluation and diagonalization of the matrix $\hat{K}$
is performed with standard numerical methods.

As the first application of ITNUMM, we solved (i) the 1D harmonic oscillator
\cite{book_Sakurai} $\mathcal{H}_{1}(q)=m\omega^{2}q^{2}/2$, using natural
units for energy and position (defined by $\hbar\omega$ and $\sqrt
{\hbar/m\omega}$, respectively), and (ii) two particles of equal mass $m$
interacting via a LJ potential $\mathcal{H}_{1}(q)=V_{\mathrm{LJ}}%
(q)\equiv\epsilon\lbrack(r_{e}/q)^{12}-2(r_{e}/q)^{6}]$. For the latter, we
used a de Boer quantum delocalization length \cite{Deckman:2009} of
$\Lambda=2^{1/6}\hbar/(r_{e}\sqrt{m\epsilon})=0.16$, corresponding to
hypothetical particles with properties between para-hydrogen ---where quantum
effects dominate---and neon---where quantum effects are present but classical
behavior dominates. In the Supplementary Material (SM), we show the grid
points, eigenvalues, and several eigenstates obtained with ITNUMM in both
cases---we also include a notebook executable in the Wolfram Research's
Mathematica software, where the interested reader can explore the technical
details of the method.

As expected, we observed that the imaginary time $\Delta\tau$ must be small
enough to reduce the relative error $\sigma$ introduced by the splitting of
the propagator---which is $\sigma\sim O(\Delta\tau^{3})$ since the second
order term vanishes for stationary states---but large enough to avoid reducing
the Gaussian elements of the $\hat{K}$ matrix to delta functions and
eventually obtaining a diagonal matrix. The latter condition is ensured by
requiring $1\gg\hat{K}_{jj}=Z_{\mathcal{H}_{1}}\left(  m/2\pi\hbar\Delta
\tau\right)  ^{D/2}/N$, which imposes a lower bound on $\Delta\tau$ for a
given $N$. In the SM, we explore the dependence of the relative error $\sigma$
on $\Delta\tau$ for a given number $N$ of grid points, and also the dependence
of $\sigma$ on $N$ in the harmonic oscillator. Remarkably, the relative error
can be fitted to $\sigma(N)\simeq0.18N^{-1.9}$, indicating a significantly
faster convergence rate than the rate expected for a MC scheme [$\sigma(N)\sim
N^{-1/2}$] \cite{McMahon:2012,Pollet:2012}. Regarding the excited states, we
found that the error becomes large for states with the highest eigenenergies.
Indeed, the number of grid points $N$ becomes insufficient to reproduce the
characteristic high frequency oscillations of wavefunctions describing highly
excited states. Yet, the agreement with exact results is very good for the
first $150$ states using $N=500$ grid points, as shown in
Fig.~\ref{fig:eneden} for the first 50 states (the whole spectrum is shown in
the SM).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{figure}
[tb]\includegraphics[width=\hsize,trim= 32 27 20 20,clip=true]{eneden3.pdf}

\caption{ (Color online) Left: Energy spectrum for $D$ LJ particles in a 1D
harmonic trap obtained with our method (circles). The exact results for $D=1$
and $2$ are shown as red solid lines. Energies are shifted to the minimum
of the potential $\min(\mathcal{H}_1)=-\protect\epsilon(D-1)D/2$. Right:
One-body densities (normalized to the number of particles) of the ground state for
distinguishable (colored lines) and indistinguishable particles (black
lines).} \label{fig:eneden}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


As a more stringent test, we now apply the method to $D$ LJ particles in a
one-dimensional harmonic trap. Potentials $V_{1}$ and $V_{2}$ are defined by
\begin{align}
V_{1}({\mathbf{q}})  &  =\sum_{\lambda=1}^{D}\frac{1}{2}m\omega^{2}q_{\lambda
}^{2},\\
V_{2}({\mathbf{q}})  &  =\sum_{\lambda<\mu}^{D}V_{\mathrm{LJ}}(|q_{\lambda
}-q_{\mu}|),
\end{align}
the de Boer length has the same value as in the example above, and $\omega
r_{e}\sqrt{m/\epsilon}=1/2$. The problem is separable only for $D=1$ or $2$,
and so a multidimensional numerical method is mandatory for $D\geq3$. In order
to reduce the number of grid points in the numerical calculation, we first
solve the problem for distinguishable particles and construct \textit{a
posteriori} the eigenstates of indistinguishable particles by symmetrizing or
anti-symmetrizing the wavefunction for spinless bosons or fermions,
respectively. Thanks to the repulsive nature of the LJ potential at short
distances we only need to evaluate $\hat{K}$ in the subspace defined by
$q_{1}>q_{2}+a,\dots,$ $q_{D-1}>q_{D}+a$, where $a$ is the core radius of the
LJ\ potential, within which the wavefunction is expected to be zero within
numerical accuracy ($a=0.63r_{e}$ in our calculations). The grid points are
sampled from the classical Boltzmann distribution of the harmonic trap in this
subspace, $p({\mathbf{q}})=Z_{V_{1}}^{-1}e^{-\Delta\tau m\omega^{2}%
|{\mathbf{q}}|^{2}/2\hbar}$ with $Z_{V_{1}}=(2\pi/\Delta\tau m\omega
^{2})^{D/2}/C_{D}(a)$, where the normalization constant obeys $C_{D}(0)=D!$.
All the two-body interactions, contained in $V_{2}({\mathbf{q}})$, are
evaluated in the matrix elements of $\hat{K}$. As mentioned above, only the
low-lying eigenstates are accurate, so we have used the Arnoldi algorithm
\cite{Arnoldi:1951} to obtain the first 50 eigenstates. We have taken into
account that many of the matrix elements are close to zero by using standard
computational techniques for sparse matrices: instead of storing the $N\times
N$ values of the matrix, only elements larger than a certain threshold were
stored. Parameters used in calculations with varying $D$ were
\begin{equation}%
\begin{array}
[c]{r|ccccc}%
D & 1 & 2 & 3 & 4 & 5\\\hline
N & 500 & 6709 & 14~394 & 36~517 & 84~690\\
\Delta\tau & 0.0055 & 0.15 & 1.5 & 1.5 & 1.5
\end{array}
\nonumber
\end{equation}
Note the relatively low total number of grid points needed to obtain results
with reasonable accuracy (a relative error of 0.002 for the $D=2$ case).
Figure~\ref{fig:states2D} shows the ground and 19th states for $D=2$ and for
the three statistics: distinguishable particles (in the above mentioned
subspace), bosons, and fermions (in the full space). The spectrum of
$\mathcal{H}$ as a function of $D$ is shown in Fig.~\ref{fig:eneden} (left
panel). We find the same spectrum for the three cases, which is a consequence
of the fermionization \cite{Paredes:2004} mechanism due to the repulsive
behavior of the LJ potential at short distances. Indeed, the bosonic and
fermionic systems show the same one-body densities in position space, as shown
in the right panel of Fig.~\ref{fig:eneden}. In all three cases the densities
show a well-defined structure, forming a quantum crystal. The displayed
one-body densities, defined as \cite{book_Lipparini}
\begin{equation}
\rho_{n}(q_{\lambda})=\int|\varphi_{n}({\mathbf{q}})|^{2}\prod_{\mu\neq
\lambda}^{D}dq_{\mu},
\end{equation}
were obtained from the nonuniform mesh as follows: first,we computed its
Fourier transform in a regular equidistant grid in momentum ($k$) space as
\begin{align}
\tilde{\rho}_{n}(k)  &  =\int|\varphi_{n}({\mathbf{q}})|^{2}e^{-ikq_{\lambda}%
}\prod_{\mu=1}^{D}dq_{\mu}\\
&  \approx\frac{Z}{N}\sum_{j=1}^{N}e^{\Delta\tau V_{1}({\mathbf{q}}_{j}%
)/\hbar-ikq_{\lambda}}|\varphi_{n}({\mathbf{q}}_{j})|^{2},
\end{align}
and then Fourier-transformed $\tilde{\rho}_{n}(k)$ back to $q_{\lambda}$-space
using standard numerical methods.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{figure}
[tb]\includegraphics[width=0.9\hsize]{states2Db.pdf}

\caption{ (Color online) Wavefunctions $\protect\varphi_n(q_1,q_2)$ of the ground and 19th states  for
$D=2$ LJ particles in a 1D harmonic trap (see text for details).
Lighter (darker) color indicates positive (negative) values of the wavefunction.
Left: distinguishable particles in the subspace $q_1>q_2$; center:
indistinguishable bosons; right: indistinguishable fermions. }
\label{fig:states2D}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


The 50 states obtained in the course of the diagonalization are sufficient to
study the behavior of the system at finite temperatures. The (unnormalized)
probability distribution of the system $p_{\beta}({\mathbf{q}})$ at finite
inverse temperature $\beta$ is defined as the thermal average
\begin{equation}
p_{\beta}({\mathbf{q}})=\sum_{n=1}^{\infty}e^{-\beta E_{n}}|\varphi
_{n}({\mathbf{q}})|^{2},
\end{equation}
and the corresponding one-body density $\rho_{\beta}(q)$ is obtained similarly
as for pure states. Figure~\ref{fig:temp} shows the one-body density for $D=4$
at three different temperatures, where the lack of structure at the highest
temperature can be understood as the melting of the quantum crystal.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{figure}
[tb]\includegraphics[width=0.9\hsize]{thermic4D.pdf}

\caption{ (Color online) One-body densities (normalized to the number of particles) for $D=4$
LJ particles in a 1D harmonic trap at three different temperatures: $\hbar
\protect\omega \protect\beta =2.8$ (dotted line), $0.7$ (dashed line), and $0.4$ (solid line).
The crystal structure disappears with increasing temperature, resulting in an
unstructured total density as in a fluid. } \label{fig:temp}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


To summarize, we have presented compelling evidence that the proposed method
achieves the original goals. Indeed, (i) the only approximation used is the
numerical discretization of space and time; (ii) the ITNUMM only requires
standard methods for sampling from an arbitrary probability distribution and
for diagonalizing real-valued sparse matrices; (iii) both ground and excited
states are obtained in the course of the diagonalization; and (iv) due to the
nonuniform nature of the grid that uses the potential to guide the sampling,
the complexity of the algorithm is significantly reduced in high-dimensional
systems. In particular, all our calculations were performed on a single
workstation with a 64-bit 2.4 GHz Quad-Core Intel Xeon E5 processor and 12 GB
of memory. Yet, the algorithm can be easily accelerated by parallelization.
The accuracy of ITNUMM can be increased by using tailored grids, larger $N$
values, or splitting methods of a higher order than in Eq. (\ref{eq:1}). In
addition to computing thermal averages---as shown here---the large set of
excited states can be also used for solving real-time quantum dynamics in a
straightforward fashion. As we have not found any \textit{a priori} limitation
to the applicability of the method, other systems described by the MDSE will
be studied in the future.

\textit{Acknowledgments.} The authors thank E. Zambrano, M. Wehrle, M. \v{S}ulc,
and F. Mazzanti for discussions. This research was supported by the Swiss NSF
NCCR MUST (Molecular Ultrafast Science \& Technology) and by the EPFL.

\bibliographystyle{apsrev4-1}
\bibliography{method_Alb}



\end{document}
tise_method_v11.bbl
%merlin.mbs apsrev4-1.bst 2010-07-25 4.21a (PWD, AO, DPC) hacked
%Control: key (0)
%Control: author (72) initials jnrlst
%Control: editor formatted (1) identically to author
%Control: production of article title (-1) disabled
%Control: page (0) single
%Control: year (1) truncated
%Control: production of eprint (0) enabled
\begin{thebibliography}{37}%
\makeatletter
\providecommand \@ifxundefined [1]{%
 \@ifx{#1\undefined}
}%
\providecommand \@ifnum [1]{%
 \ifnum #1\expandafter \@firstoftwo
 \else \expandafter \@secondoftwo
 \fi
}%
\providecommand \@ifx [1]{%
 \ifx #1\expandafter \@firstoftwo
 \else \expandafter \@secondoftwo
 \fi
}%
\providecommand \natexlab [1]{#1}%
\providecommand \enquote  [1]{``#1''}%
\providecommand \bibnamefont  [1]{#1}%
\providecommand \bibfnamefont [1]{#1}%
\providecommand \citenamefont [1]{#1}%
\providecommand \href@noop [0]{\@secondoftwo}%
\providecommand \href [0]{\begingroup \@sanitize@url \@href}%
\providecommand \@href[1]{\@@startlink{#1}\@@href}%
\providecommand \@@href[1]{\endgroup#1\@@endlink}%
\providecommand \@sanitize@url [0]{\catcode `\\12\catcode `\$12\catcode
  `\&12\catcode `\#12\catcode `\^12\catcode `\_12\catcode `\%12\relax}%
\providecommand \@@startlink[1]{}%
\providecommand \@@endlink[0]{}%
\providecommand \url  [0]{\begingroup\@sanitize@url \@url }%
\providecommand \@url [1]{\endgroup\@href {#1}{\urlprefix }}%
\providecommand \urlprefix  [0]{URL }%
\providecommand \Eprint [0]{\href }%
\providecommand \doibase [0]{http://dx.doi.org/}%
\providecommand \selectlanguage [0]{\@gobble}%
\providecommand \bibinfo  [0]{\@secondoftwo}%
\providecommand \bibfield  [0]{\@secondoftwo}%
\providecommand \translation [1]{[#1]}%
\providecommand \BibitemOpen [0]{}%
\providecommand \bibitemStop [0]{}%
\providecommand \bibitemNoStop [0]{.\EOS\space}%
\providecommand \EOS [0]{\spacefactor3000\relax}%
\providecommand \BibitemShut  [1]{\csname bibitem#1\endcsname}%
\let\auto@bib@innerbib\@empty
%</preamble>
\bibitem [{\citenamefont {Harris}\ \emph {et~al.}(1965)\citenamefont {Harris},
  \citenamefont {Engerholm},\ and\ \citenamefont {Gwinn}}]{Harris:1965}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~O.}\ \bibnamefont
  {Harris}}, \bibinfo {author} {\bibfnamefont {G.~G.}\ \bibnamefont
  {Engerholm}}, \ and\ \bibinfo {author} {\bibfnamefont {W.~D.}\ \bibnamefont
  {Gwinn}},\ }\href {\doibase 110.1063/1.1696963} {\bibfield  {journal}
  {\bibinfo  {journal} {J. Chem. Phys.}\ }\textbf {\bibinfo {volume} {43}},\
  \bibinfo {pages} {1515} (\bibinfo {year} {1965})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kosloff}(1988)}]{Kosloff:1988}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
  {Kosloff}},\ }\href {\doibase 10.1021/j100319a003} {\bibfield  {journal}
  {\bibinfo  {journal} {J. Phys. Chem.}\ }\textbf {\bibinfo {volume} {92}},\
  \bibinfo {pages} {2087} (\bibinfo {year} {1988})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Shimshovitz}\ and\ \citenamefont
  {Tannor}(2012)}]{Shimshovitz:2012}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Shimshovitz}}\ and\ \bibinfo {author} {\bibfnamefont {D.~J.}\ \bibnamefont
  {Tannor}},\ }\href {\doibase 10.1103/PhysRevLett.109.070402} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo
  {volume} {109}},\ \bibinfo {pages} {070402} (\bibinfo {year}
  {2012})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Pollet}(2012)}]{Pollet:2012}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {L.}~\bibnamefont
  {Pollet}},\ }\href {\doibase 10.1088/0034-4885/75/9/094501} {\bibfield
  {journal} {\bibinfo  {journal} {Rep. Prog. Phys.}\ }\textbf {\bibinfo
  {volume} {75}},\ \bibinfo {pages} {094501} (\bibinfo {year}
  {2012})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Parr}\ and\ \citenamefont
  {Yang}(1994)}]{book_Parr_Yang}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.~G.}\ \bibnamefont
  {Parr}}\ and\ \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont {Yang}},\
  }\href@noop {} {\emph {\bibinfo {title} {{Density-Functional Theory of Atoms
  and Molecules}}}},\ \bibinfo {edition} {{2}}\ ed.\ (\bibinfo  {publisher}
  {{Oxford University Press}},\ \bibinfo {year} {1994})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Dalfovo}\ \emph {et~al.}(1995)\citenamefont
  {Dalfovo}, \citenamefont {Lastri}, \citenamefont {Pricaupenko}, \citenamefont
  {Stringari},\ and\ \citenamefont {Treiner}}]{Dalfovo:1995}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {F.}~\bibnamefont
  {Dalfovo}}, \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Lastri}},
  \bibinfo {author} {\bibfnamefont {L.}~\bibnamefont {Pricaupenko}}, \bibinfo
  {author} {\bibfnamefont {S.}~\bibnamefont {Stringari}}, \ and\ \bibinfo
  {author} {\bibfnamefont {J.}~\bibnamefont {Treiner}},\ }\href {\doibase
  10.1103/PhysRevB.52.1193} {\bibfield  {journal} {\bibinfo  {journal} {Phys.
  Rev. B}\ }\textbf {\bibinfo {volume} {52}},\ \bibinfo {pages} {1193}
  (\bibinfo {year} {1995})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Savenko}\ \emph {et~al.}(2013)\citenamefont
  {Savenko}, \citenamefont {Liew},\ and\ \citenamefont
  {Shelykh}}]{Savenko:2013}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {I.~G.}\ \bibnamefont
  {Savenko}}, \bibinfo {author} {\bibfnamefont {T.~C.~H.}\ \bibnamefont
  {Liew}}, \ and\ \bibinfo {author} {\bibfnamefont {I.~A.}\ \bibnamefont
  {Shelykh}},\ }\href {\doibase 10.1103/PhysRevLett.110.127402} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo
  {volume} {110}},\ \bibinfo {pages} {127402} (\bibinfo {year}
  {2013})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Rogel-Salazar}(2013)}]{Rogel-Salazar:2013}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Rogel-Salazar}},\ }\href {\doibase 10.1088/0143-0807/34/2/247} {\bibfield
  {journal} {\bibinfo  {journal} {Eur. J. Phys.}\ }\textbf {\bibinfo {volume}
  {34}},\ \bibinfo {pages} {247} (\bibinfo {year} {2013})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Billing}\ and\ \citenamefont
  {Adhikari}(2000)}]{Billing:2000}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Billing}}\ and\ \bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Adhikari}},\ }\href {\doibase 10.1016/S0009-2614(00)00358-4} {\bibfield
  {journal} {\bibinfo  {journal} {Chem. Phys. Lett.}\ }\textbf {\bibinfo
  {volume} {321}},\ \bibinfo {pages} {197} (\bibinfo {year}
  {2000})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Runge}\ and\ \citenamefont
  {Gross}(1984)}]{Runge_Gross:1984}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Runge}}\ and\ \bibinfo {author} {\bibfnamefont {E.~K.~U.}\ \bibnamefont
  {Gross}},\ }\href {\doibase 10.1103/PhysRevLett.52.997} {\bibfield  {journal}
  {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo {volume} {52}},\
  \bibinfo {pages} {997} (\bibinfo {year} {1984})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Abid}\ \emph {et~al.}(2003)\citenamefont {Abid},
  \citenamefont {Huepe}, \citenamefont {Metens}, \citenamefont {Nore},
  \citenamefont {Pham}, \citenamefont {Tuckerman},\ and\ \citenamefont
  {Brachet}}]{Abid:2003}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Abid}}, \bibinfo {author} {\bibfnamefont {C.}~\bibnamefont {Huepe}},
  \bibinfo {author} {\bibfnamefont {S.}~\bibnamefont {Metens}}, \bibinfo
  {author} {\bibfnamefont {C.}~\bibnamefont {Nore}}, \bibinfo {author}
  {\bibfnamefont {C.~T.}\ \bibnamefont {Pham}}, \bibinfo {author}
  {\bibfnamefont {L.~S.}\ \bibnamefont {Tuckerman}}, \ and\ \bibinfo {author}
  {\bibfnamefont {M.~E.}\ \bibnamefont {Brachet}},\ }\href {\doibase
  10.1016/j.fluiddyn.2003.09.001} {\bibfield  {journal} {\bibinfo  {journal}
  {Fluid Dyn. Res.}\ }\textbf {\bibinfo {volume} {33}},\ \bibinfo {pages} {509}
  (\bibinfo {year} {2003})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Wyatt}\ and\ \citenamefont
  {Trahan}(2005)}]{book_Wyatt}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.~E.}\ \bibnamefont
  {Wyatt}}\ and\ \bibinfo {author} {\bibfnamefont {C.~J.}\ \bibnamefont
  {Trahan}},\ }\href@noop {} {\emph {\bibinfo {title} {{Quantum Dynamics with
  Trajectories: Introduction to Quantum Hydrodynamics}}}},\ \bibinfo {edition}
  {1st}\ ed.\ (\bibinfo  {publisher} {{Addison-Wesley}},\ \bibinfo {year}
  {2005})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Alvarez}\ \emph {et~al.}(2011)\citenamefont
  {Alvarez}, \citenamefont {da~Silva}, \citenamefont {Ponce},\ and\
  \citenamefont {Dagotto}}]{Alvarez:2011}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Alvarez}}, \bibinfo {author} {\bibfnamefont {L.~G. G. V.~D.}\ \bibnamefont
  {da~Silva}}, \bibinfo {author} {\bibfnamefont {E.}~\bibnamefont {Ponce}}, \
  and\ \bibinfo {author} {\bibfnamefont {E.}~\bibnamefont {Dagotto}},\ }\href
  {\doibase 10.1103/PhysRevE.84.056706} {\bibfield  {journal} {\bibinfo
  {journal} {Phys. Rev. E}\ }\textbf {\bibinfo {volume} {84}},\ \bibinfo
  {pages} {056706} (\bibinfo {year} {2011})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Tkatchenko}\ \emph {et~al.}(2012)\citenamefont
  {Tkatchenko}, \citenamefont {DiStasio}, \citenamefont {Car},\ and\
  \citenamefont {Scheffler}}]{Tkatchenko:2012}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Tkatchenko}}, \bibinfo {author} {\bibfnamefont {R.~A.}\ \bibnamefont
  {DiStasio}}, \bibinfo {author} {\bibfnamefont {R.}~\bibnamefont {Car}}, \
  and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Scheffler}},\ }\href
  {http://link.aps.org/doi/10.1103/PhysRevLett.108.236402} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo
  {volume} {108}},\ \bibinfo {pages} {236402} (\bibinfo {year}
  {2012})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Mart\'in-Delgado}\ and\ \citenamefont
  {Sierra}(1996)}]{Delgado:1996}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.~A.}\ \bibnamefont
  {Mart\'in-Delgado}}\ and\ \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Sierra}},\ }\href {http://link.aps.org/doi/10.1103/PhysRevLett.76.1146}
  {\bibfield  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf
  {\bibinfo {volume} {76}},\ \bibinfo {pages} {1146} (\bibinfo {year}
  {1996})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Corney}\ and\ \citenamefont
  {Drummond}(2004)}]{Corney:2004}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~F.}\ \bibnamefont
  {Corney}}\ and\ \bibinfo {author} {\bibfnamefont {P.~D.}\ \bibnamefont
  {Drummond}},\ }\href {http://link.aps.org/doi/10.1103/PhysRevLett.93.260401}
  {\bibfield  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf
  {\bibinfo {volume} {93}},\ \bibinfo {pages} {260401} (\bibinfo {year}
  {2004})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Torres-Vega}\ and\ \citenamefont
  {Frederick}(1991)}]{Torres-Vega:1991}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Torres-Vega}}\ and\ \bibinfo {author} {\bibfnamefont {J.~H.}\ \bibnamefont
  {Frederick}},\ }\href {http://link.aps.org/doi/10.1103/PhysRevLett.67.2601}
  {\bibfield  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf
  {\bibinfo {volume} {67}},\ \bibinfo {pages} {2601} (\bibinfo {year}
  {1991})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {McMahon}\ \emph {et~al.}(2012)\citenamefont
  {McMahon}, \citenamefont {Morales}, \citenamefont {Pierleoni},\ and\
  \citenamefont {Ceperley}}]{McMahon:2012}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~M.}\ \bibnamefont
  {McMahon}}, \bibinfo {author} {\bibfnamefont {M.~A.}\ \bibnamefont
  {Morales}}, \bibinfo {author} {\bibfnamefont {C.}~\bibnamefont {Pierleoni}},
  \ and\ \bibinfo {author} {\bibfnamefont {D.~M.}\ \bibnamefont {Ceperley}},\
  }\href {\doibase 10.1103/RevModPhys.84.1607} {\bibfield  {journal} {\bibinfo
  {journal} {Rev. Mod. Phys.}\ }\textbf {\bibinfo {volume} {84}},\ \bibinfo
  {pages} {1607} (\bibinfo {year} {2012})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kapitza}(1938)}]{Kapitza:1938}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
  {Kapitza}},\ }\href {\doibase 10.1038/141074a0} {\bibfield  {journal}
  {\bibinfo  {journal} {Nature}\ }\textbf {\bibinfo {volume} {141}},\ \bibinfo
  {pages} {74} (\bibinfo {year} {1938})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Leggett}(1972)}]{Leggett:1972}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.~J.}\ \bibnamefont
  {Leggett}},\ }\href {\doibase 10.1103/PhysRevLett.29.1227} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo
  {volume} {29}},\ \bibinfo {pages} {1227} (\bibinfo {year}
  {1972})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Grebenev}\ \emph {et~al.}(2000)\citenamefont
  {Grebenev}, \citenamefont {Sartakov}, \citenamefont {Toennies},\ and\
  \citenamefont {Vilesov}}]{Grebenev:2000}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Grebenev}}, \bibinfo {author} {\bibfnamefont {B.}~\bibnamefont {Sartakov}},
  \bibinfo {author} {\bibfnamefont {J.~P.}\ \bibnamefont {Toennies}}, \ and\
  \bibinfo {author} {\bibfnamefont {A.~F.}\ \bibnamefont {Vilesov}},\ }\href
  {\doibase 10.1126/science.289.5484.1532} {\bibfield  {journal} {\bibinfo
  {journal} {Science}\ }\textbf {\bibinfo {volume} {289}},\ \bibinfo {pages}
  {1532} (\bibinfo {year} {2000})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Anglin}\ and\ \citenamefont
  {Ketterle}(2002)}]{Anglin_Ketterle:2002}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~R.}\ \bibnamefont
  {Anglin}}\ and\ \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont
  {Ketterle}},\ }\href {\doibase 10.1038/416211a} {\bibfield  {journal}
  {\bibinfo  {journal} {Nature}\ }\textbf {\bibinfo {volume} {416}},\ \bibinfo
  {pages} {211} (\bibinfo {year} {2002})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Balibar}(2010)}]{Balibar:2010}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Balibar}},\ }\href {\doibase 10.1038/nature08913} {\bibfield  {journal}
  {\bibinfo  {journal} {Nature}\ }\textbf {\bibinfo {volume} {464}},\ \bibinfo
  {pages} {176} (\bibinfo {year} {2010})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Bloch}(2005)}]{Bloch:2005}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {I.}~\bibnamefont
  {Bloch}},\ }\href {\doibase 10.1038/nphys138} {\bibfield  {journal} {\bibinfo
   {journal} {Nature Physics}\ }\textbf {\bibinfo {volume} {1}},\ \bibinfo
  {pages} {23} (\bibinfo {year} {2005})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ronzheimer}\ \emph {et~al.}(2013)\citenamefont
  {Ronzheimer}, \citenamefont {Schreiber}, \citenamefont {Braun}, \citenamefont
  {Hodgman}, \citenamefont {Langer}, \citenamefont {McCulloch}, \citenamefont
  {Heidrich-Meisner}, \citenamefont {Bloch},\ and\ \citenamefont
  {Schneider}}]{Ronzheimer:2013}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~P.}\ \bibnamefont
  {Ronzheimer}}, \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Schreiber}}, \bibinfo {author} {\bibfnamefont {S.}~\bibnamefont {Braun}},
  \bibinfo {author} {\bibfnamefont {S.~S.}\ \bibnamefont {Hodgman}}, \bibinfo
  {author} {\bibfnamefont {S.}~\bibnamefont {Langer}}, \bibinfo {author}
  {\bibfnamefont {I.~P.}\ \bibnamefont {McCulloch}}, \bibinfo {author}
  {\bibfnamefont {F.}~\bibnamefont {Heidrich-Meisner}}, \bibinfo {author}
  {\bibfnamefont {I.}~\bibnamefont {Bloch}}, \ and\ \bibinfo {author}
  {\bibfnamefont {U.}~\bibnamefont {Schneider}},\ }\href
  {http://link.aps.org/doi/10.1103/PhysRevLett.110.205301} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo
  {volume} {110}},\ \bibinfo {pages} {205301} (\bibinfo {year}
  {2013})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Panfil}\ \emph {et~al.}(2013)\citenamefont {Panfil},
  \citenamefont {Nardis},\ and\ \citenamefont {Caux}}]{Panfil:2013}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Panfil}}, \bibinfo {author} {\bibfnamefont {J.~D.}\ \bibnamefont {Nardis}},
  \ and\ \bibinfo {author} {\bibfnamefont {J.-S.}\ \bibnamefont {Caux}},\
  }\href {http://link.aps.org/doi/10.1103/PhysRevLett.110.125302} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo
  {volume} {110}},\ \bibinfo {pages} {125302} (\bibinfo {year}
  {2013})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Vignolo}\ and\ \citenamefont
  {Minguzzi}(2013)}]{Vignolo:2013}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
  {Vignolo}}\ and\ \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Minguzzi}},\ }\href {http://link.aps.org/doi/10.1103/PhysRevLett.110.020403}
  {\bibfield  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf
  {\bibinfo {volume} {110}},\ \bibinfo {pages} {020403} (\bibinfo {year}
  {2013})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Gring}\ \emph {et~al.}(2012)\citenamefont {Gring},
  \citenamefont {Kuhnert}, \citenamefont {Langen}, \citenamefont {Kitagawa},
  \citenamefont {Rauer}, \citenamefont {Schreitl}, \citenamefont {Mazets},
  \citenamefont {Smith}, \citenamefont {Demler},\ and\ \citenamefont
  {Schmiedmayer}}]{Gring:2012}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Gring}}, \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Kuhnert}},
  \bibinfo {author} {\bibfnamefont {T.}~\bibnamefont {Langen}}, \bibinfo
  {author} {\bibfnamefont {T.}~\bibnamefont {Kitagawa}}, \bibinfo {author}
  {\bibfnamefont {B.}~\bibnamefont {Rauer}}, \bibinfo {author} {\bibfnamefont
  {M.}~\bibnamefont {Schreitl}}, \bibinfo {author} {\bibfnamefont
  {I.}~\bibnamefont {Mazets}}, \bibinfo {author} {\bibfnamefont {D.~A.}\
  \bibnamefont {Smith}}, \bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Demler}}, \ and\ \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Schmiedmayer}},\ }\href {\doibase 10.1126/science.1224953} {\bibfield
  {journal} {\bibinfo  {journal} {Science}\ }\textbf {\bibinfo {volume}
  {337}},\ \bibinfo {pages} {1318} (\bibinfo {year} {2012})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Kinoshita}\ \emph {et~al.}(2006)\citenamefont
  {Kinoshita}, \citenamefont {Wenger},\ and\ \citenamefont
  {Weiss}}]{Kinoshita:2006}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
  {Kinoshita}}, \bibinfo {author} {\bibfnamefont {T.}~\bibnamefont {Wenger}}, \
  and\ \bibinfo {author} {\bibfnamefont {D.~S.}\ \bibnamefont {Weiss}},\ }\href
  {\doibase 10.1038/nature04693} {\bibfield  {journal} {\bibinfo  {journal}
  {Nature}\ }\textbf {\bibinfo {volume} {440}},\ \bibinfo {pages} {900}
  (\bibinfo {year} {2006})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Paredes}\ \emph {et~al.}(2004)\citenamefont
  {Paredes}, \citenamefont {Widera}, \citenamefont {Murg}, \citenamefont
  {Mandel}, \citenamefont {Folling}, \citenamefont {Cirac}, \citenamefont
  {Shlyapnikov}, \citenamefont {Hansch},\ and\ \citenamefont
  {Bloch}}]{Paredes:2004}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {B.}~\bibnamefont
  {Paredes}}, \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Widera}},
  \bibinfo {author} {\bibfnamefont {V.}~\bibnamefont {Murg}}, \bibinfo {author}
  {\bibfnamefont {O.}~\bibnamefont {Mandel}}, \bibinfo {author} {\bibfnamefont
  {S.}~\bibnamefont {Folling}}, \bibinfo {author} {\bibfnamefont
  {I.}~\bibnamefont {Cirac}}, \bibinfo {author} {\bibfnamefont {G.~V.}\
  \bibnamefont {Shlyapnikov}}, \bibinfo {author} {\bibfnamefont {T.~W.}\
  \bibnamefont {Hansch}}, \ and\ \bibinfo {author} {\bibfnamefont
  {I.}~\bibnamefont {Bloch}},\ }\href {\doibase 10.1038/nature02530} {\bibfield
   {journal} {\bibinfo  {journal} {Nature}\ }\textbf {\bibinfo {volume}
  {429}},\ \bibinfo {pages} {277} (\bibinfo {year} {2004})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Bloch}\ \emph {et~al.}(2008)\citenamefont {Bloch},
  \citenamefont {Dalibard},\ and\ \citenamefont {Zwerger}}]{Bloch:2008}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {I.}~\bibnamefont
  {Bloch}}, \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {Dalibard}}, \
  and\ \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont {Zwerger}},\ }\href
  {\doibase 10.1103/RevModPhys.80.885} {\bibfield  {journal} {\bibinfo
  {journal} {Rev. Mod. Phys.}\ }\textbf {\bibinfo {volume} {80}},\ \bibinfo
  {pages} {885} (\bibinfo {year} {2008})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Morton}\ and\ \citenamefont
  {Mayers}(2005)}]{book_Morton_Mayers}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.}~\bibnamefont
  {Morton}}\ and\ \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Mayers}},\ }\href@noop {} {\emph {\bibinfo {title} {{Numerical Solution of
  Partial Differential Equations, An Introduction}}}},\ \bibinfo {edition}
  {{2}}\ ed.\ (\bibinfo  {publisher} {{Cambridge University Press}},\ \bibinfo
  {year} {2005})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Popov}(2005)}]{Popov:2005}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {V.~S.}\ \bibnamefont
  {Popov}},\ }\href {\doibase 10.1134/1.1903097} {\bibfield  {journal}
  {\bibinfo  {journal} {Phys. Atom. Nuclei}\ }\textbf {\bibinfo {volume}
  {68}},\ \bibinfo {pages} {686} (\bibinfo {year} {2005})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Sakurai}(1993)}]{book_Sakurai}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~J.}\ \bibnamefont
  {Sakurai}},\ }\href@noop {} {\emph {\bibinfo {title} {{Modern Quantum
  Mechanics }}}},\ \bibinfo {edition} {{1}}\ ed.\ (\bibinfo  {publisher}
  {{Addison-Wesley}},\ \bibinfo {year} {1993})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Deckman}\ and\ \citenamefont
  {Mandelshtam}(2009)}]{Deckman:2009}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Deckman}}\ and\ \bibinfo {author} {\bibfnamefont {V.~A.}\ \bibnamefont
  {Mandelshtam}},\ }\href {\doibase 10.1021/jp900095f} {\bibfield  {journal}
  {\bibinfo  {journal} {J. Phys. Chem. A}\ }\textbf {\bibinfo {volume} {26}},\
  \bibinfo {pages} {7394} (\bibinfo {year} {2009})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Arnoldi}(1951)}]{Arnoldi:1951}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.~E.}\ \bibnamefont
  {Arnoldi}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Q.
  Appl. Math}\ }\textbf {\bibinfo {volume} {9}},\ \bibinfo {pages} {17}
  (\bibinfo {year} {1951})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Lipparini}(2008)}]{book_Lipparini}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Lipparini}},\ }\href@noop {} {\emph {\bibinfo {title} {{Modern Many-Particle
  Physics: Atomic Gases, Nanostructures and Quantum Liquids}}}}\ (\bibinfo
  {publisher} {{World Scientific}},\ \bibinfo {year} {2008})\BibitemShut
  {NoStop}%
\end{thebibliography}%
