# Signatures of the BCS-BEC crossover in the yrast spectra of Fermi quantum
rings

Ulrich Ebling<EMAIL_ADDRESS>Centre for Theoretical Chemistry and
Physics, New Zealand Institute for Advanced Study, Massey University, Private
Bag 902104, North Shore, Auckland 0745, New Zealand  Dodd-Walls Centre for
Photonic and Quantum Technologies, PO Box 56, Dunedin 9056, New Zealand  Ali
Alavi  Max Planck Institute for Solid State Research, Heisenbergstraße 1,
70569 Stuttgart, Germany  Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge, CB2 1EW, United Kingdom  Joachim Brand  Centre for
Theoretical Chemistry and Physics, New Zealand Institute for Advanced Study,
Massey University, Private Bag 902104, North Shore, Auckland 0745, New Zealand
Dodd-Walls Centre for Photonic and Quantum Technologies, PO Box 56, Dunedin
9056, New Zealand  Max Planck Institute for Solid State Research,
Heisenbergstraße 1, 70569 Stuttgart, Germany

###### Abstract

We study properties of the lowest energy states at non-zero total momentum
(yrast states) of the Hubbard model for spin-$\frac{1}{2}$ fermions in the
quantum ring configuration with attractive on-site interaction at low density.
In the one-dimensional (1D) case we solve the Hubbard model using the Bethe
ansatz, while for the crossover into the 2D regime we use the Full-
Configuration-Interaction Quantum Monte-Carlo method (FCIQMC) to obtain the
yrast states for the spin-balanced Fermi system. We show how the yrast
excitation spectrum changes from the 1D to the 2D regime and how pairing
affects the yrast spectra. We also find signatures of fragmented condensation
for certain yrast states usually associated with dark solitons.

##  I Introduction

The crossover from a fermionic superfluid of weakly-bound Cooper pairs (BCS
regime) to a Bose-Einstein condensate (BEC) of strongly-bound dimers is a
paradigmatic quantum many-body problem [1, 2, 3, 4]. Our understanding of this
problem is still limited, as strong quantum correlations and the absence of a
small parameter pose severe challenges for theoretical and computational
approaches. While bulk systems have been studied extensively in recent years
using theory [3, 5, 6] and experiments with quantum gases [7, 8, 9, 10, 11,
12, 13, 6, 14, 15], the advent of quantum gas microscopes [16, 17, 18, 19] and
micro traps [20, 21, 22, 23] has opened up the opportunity to experiment with
systems that are small enough to perform exact numerical calculations on.

Of particular interest are ring configurations, where translational invariance
along one spatial dimension makes (angular) momentum a good quantum number.
This allows for the study of yrast states, which are defined as the lowest
energy state at given value of the total momentum. Yrast states in a bosonic
superfluid are intimately connected [24, 25, 26, 27, 28, 29, 30, 31] to
localized nonlinear waves known as dark solitons [32]. It was shown that
measuring the position of all bosons in an yrast state reveals a dark-soliton-
like particle depletion [30], and that wave-packet-like superpositions of
yrast states emulate the behavior of classical dark solitons [31]. While the
yrast states are fragmented quantum condensates, breaking the translational
symmetry restores single condensation and classical soliton features [28].
Dark solitons in Fermi superfluids have been identified in experiments [33,
34, 35], but many predictions from mean-field theory have not yet been tested
[36, 37, 38, 39, 40, 41, 42]. Moreover, there is an intriguing connection [43,
44] between dark solitons and the predicted Fulde-Ferrell-Larkin-Ovchinnikov
(FFLO) phase of imbalanced superfluids [45, 46]. Dispersion relations of yrast
states were analyzed in the context of dark soliton physics in the Yang-Gaudin
model, a one-dimensional (1D) Bethe-ansatz solvable model of a fermionic
superfluid in Ref. [47]. An overview of computational studies of quantum rings
can be found in the recent review literature [48, 49, 50].

Beyond the purely one-dimensional models of quantum rings, a second spatial
dimension can be added by considering stripe- or ladder-type lattice
configurations as it is done in this work. In Refs. [48, 49] these are
referred to as quasi-1D geometries. While the Bethe ansatz is unavailable for
such models and mean-field theory is not valid in the strongly-correlated
regime, direct numerical simulation is very challenging due to the fact that
Hilbert space size increases exponentially both with particle number and the
number of lattice sites. Quantum Monte Carlo (QMC) simulations are still
possible, although the _fermion sign problem_ [51, 52] provides a challenge
for the simulation of fermionic many-body problems.

QMC methods relevant in the field of ultracold quantum gases each have their
own strengths and weaknesses. Diffusion Monte Carlo has no basis set
dependence but either converges to a bosonic ground state or requires node-
fixing, which introduces an approximation [53]. Auxiliary-field QMC is sign-
problem free for the attractive balanced Hubbard model [54], but the Hubbard-
Stratonovich transformation involved in this method breaks symmetries of the
Hamiltonian and thus does not allow for the study of yrast states. Recent work
has suggested a solution [55] but it has yet to be seen whether the method can
be implemented efficiently. Determinant Monte Carlo and related methods can
handle finite temperature and extrapolate to zero temperature, but they fix
the chemical potential instead of particle number [56, 57, 58, 59]. All these
existing methods have in common that they can study overall ground-state
properties while it is not possible to study yrast states, because the total
momentum cannot be easily constrained.

For this work we use Full-Configuration-Interaction Quantum Monte Carlo
(FCIQMC), a method originally developed for strongly correlated electrons in
the context of quantum chemistry [60, 61]. FCIQMC has been applied with great
success to a large number of problems in this field [62, 63], and recently to
ultracold atoms [64, 65]. This method can find the ground state energy and
many-body wave function in a Fermi system by expanding the wave function into
a set of Slater determinants. A stochastic version of exact diagonalization of
the Hamiltonian in this basis is achieved by simulating the dynamics of a
walker population in Slater determinant space. FCIQMC mitigates the sign
problem by walker anihilation to a certain degree but does not eliminate it
[66]. By performing a stochastic projection to the ground state of a
Hamiltonian in a given Fock basis directly, i.e. without resorting to a
Hubbard-Stratonovich transformation, it is easy to respect symmetries of the
Hamiltonian. In particular, it is possible to obtain energies and observables
from yrast states by ground state projection in a plane-wave basis because the
FCIQMC algorithm conserves total momentum if the Hamiltonian does. With the
FCI method taking into account all correlations in the system, we probe the
BCS-BEC crossover from the non-interacting to the strongly attractive regime
in the Hubbard model.

In this paper, we present a QMC study of yrast states in the Hubbard model for
spin-$\frac{1}{2}$ fermions. Using a filling factor much smaller than one,
this system resembles a continuum superfluid with the difference that momentum
is replaced by lattice momentum. We study the crossover from 1D to 2D geometry
in the case of attractive on-site interactions. For 1D Hubbard chains, we
obtain exact results using the Bethe ansatz and compare them to QMC results.
Then we use FCIQMC to investigate the crossover into the 2D geometry by
increasing the number of sites in the transverse direction. We find that in
2D, the shape of the yrast dispersion changes considerably, but that for
increasing interaction strength, the more typical shape of the 1D spectrum is
restored, which indicates soliton-like physics and is a signature of the
transition into the superfluid regime. We investigate in more detail the
behavior of the local minima of the yrast spectrum at the so-called “umklapp”
points where sufficient quasi-momentum is added to boost either all or half of
the constituent fermions by one unit in order to form a ring current. We find
signatures of the transition from non-interacting Fermi gas to a paired
superfluid and of the BCS to BEC crossover in the excitation energy and in the
pair correlation functions for the first half and full umklapp points. Last,
we focus on the yrast states around the maxima of the dispersion, which are
related to dark solitons. We calculate the inertial mass of possible solitons
for different system sizes and interaction strengths. We find an increase of
the inertial mass by a factor of 2 when the transverse dimension is large
enough for the system to be considered truly 2D, which is indicative of a
transition from dark soliton to a solitonic vortex [67, 68]. From mean-field
and basic hydrodynamic theory, in the 1D to 2D crossover dark solitons are
replaced as stable yrast excitations by solitonic vortices [68, 69], or vortex
pairs [41, 70], which have larger inertial mass [34, 71, 72]. By looking at
the pair densities of these yrast states, we find that fragmented condensation
into more than one momentum state takes place, as expected for superfluid
yrast states [28].

This paper is organized as follows: After introducing the model in Sec. II and
the FCIQMC approach in Sec. III, we discuss yrast spectra of a 1D Hubbard
chain obtained by the Bethe ansatz in Sec. IV. The energies and spin
correlation functions for the umklapp points in the 1D to 2D crossover are
discussed in Sec. V before analysing the physics of the maxima of the yrast
dispersion by computing their effective mass and momentum-space pair densities
in Sec. VI, and drawing conclusions in Sec. VII.

##  II System

To study yrast states in an ultracold fermionic superfluid, we use the Hubbard
model in the regime of low densities. In this regime, the Hubbard model
approximates a discretized free space. The correspondence to free space
becomes exact in the low-density limit. We focus on the crossover between the
1D and 2D geometry, therefore our Hubbard model corresponds to a rectangular
lattice with $L\times W$ sites, where $1\leq W<L$, but with the same lattice
spacing $\alpha$ in both dimensions. We are using periodic boundary conditions
in both directions, and thus our systems has the topology of a torus. The
Hubbard Hamiltonian in momentum representation is

$H=\sum_{\vec{k},s}\epsilon_{\vec{k}}c_{\vec{k}s}^{\dagger}c_{\vec{k}s}+\frac{U}{LW}\sum_{\vec{k}_{1},\vec{k}_{2},\vec{k}_{3}}c_{\vec{k}_{1}\uparrow}^{\dagger}c_{\vec{k}_{2}\downarrow}^{\dagger}c_{\vec{k}_{3}\downarrow}c_{\vec{k}_{1}+\vec{k}_{2}-\vec{k}_{3}\uparrow},$
(1)

where the operators $c^{(\dagger)}_{\vec{k},s}$ create (annihilate) a fermion
with lattice momentum $\hbar\vec{k}$ and spin $s$. Wave vectors can take on
values $\vec{k}=(k_{x},k_{y})=2\pi(n_{x}/\alpha L,n_{y}/\alpha W)$, where
$n_{x},n_{y}$ are integers. The single-particle dispersion for the Hubbard
model is given by

$\epsilon_{\vec{k}}=2t(2-\cos(k_{x}\alpha)-\cos(k_{y}\alpha))$ (2)

where $t$ is the hopping amplitude and $U<0$ is the interaction parameter. In
the low density regime, mostly the low-lying momentum states are occupied
where the dispersion relation (2) is approximately parabolic. Thus the Hubbard
model approximates a continuum Fermi gas.

Throughout this paper we present results obtained for a particle number of
$N=10$, with 5 fermions in each spin state, and a lattice length of $L=21$,
while the width $W$ varies from 1 to 11. Energies will be given in units of
hopping amplitude $t$ and momenta in longitudinal lattice units
$P_{0}=\frac{2\pi\hbar}{\alpha L}$.

##  III FCIQMC

FCIQMC is a numerical method originally created in the context of strongly
correlated electron systems and quantum chemistry. It can find ground states
of fermionic many-body systems by expanding the many-body wave function in
terms of Slater determinants

$|\Psi\rangle=\sum_{i}C_{i}|D_{i}\rangle,$ (3)

which in our case are of the form

$|D\rangle=c_{\vec{k}_{1},s_{1}}^{\dagger}\ldots
c_{\vec{k}_{N},s_{N}}^{\dagger}|\text{vac}\rangle.$ (4)

FCIQMC then obtains the expansion coefficients $C_{i}$ by using a stochastic
population dynamics approach to solve the imaginary-time Schrödinger equation.
The automatic antisymmetrization of the wave function by expanding it in
Slater determinants ensures that unlike Diffusion Monte-Carlo, FCIQMC always
finds a fermionic wave function. FCIQMC’s capacity for overcoming the so-
called “fermion sign problem”, which here manifests in fluctuations of the
sign of each of the coefficients $C_{i}$ and which cannot be pre-determined,
depends on the importance of annihilation events among the walkers in
establishing the sign structure of the sampled wavefunction. When this effect
is strong, the full FCIQMC method requires a number of walkers which scales
with the size of the Hilbert space, making it impractical for large spaces.

To counter this, we use a range of modifications to FCIQMC, which facilitates
calculations when the sign problem is strong. We make use of a similarity-
transformed Hamiltonian which makes the many-body wave function more compact
in Hilbert space [73]. We also use the initiator approximation [61], which can
introduce a bias into the energy that disappears in the limit of large walker
number. In order to control this undesirable bias, we first compare QMC
results with exact results in the 1D case and then adjust the walker number
until the initiator bias is eliminated. For the 2D systems, we successively
increase the walker number with increasing $W$. Walker numbers used in this
paper range from $N_{W}=1\times 10^{6}$ for weakly-interacting 1D chains to
$N_{W}=2\times 10^{8}$ for 2D systems at $U/t=-5$. The most demanding
computations were run on up to 400 processor cores using up to 2 Gigabyte
memory per core.

##  IV Bethe ansatz results

For a one-dimensional Hubbard chain ($W=1$), the system is integrable and the
Hamiltonian (1) can be diagonalized using the Bethe ansatz [74]. For a
balanced Fermi system with $N$ fermions, energy and total momentum are given
by

$\displaystyle E$ $\displaystyle=-2t\sum_{j=1}^{N}\cos(\kappa_{j})+U(L-2N),$
(5) $\displaystyle P/P_{0}$
$\displaystyle=\sum_{j=1}^{N}\kappa_{j}\quad\text{mod}\quad 2\pi,$ (6)

with $N$ dimensionless quasi-momenta $\kappa_{j}$ which must be obtained,
alongside $N/2$ rapidities $\Lambda_{\alpha}$, by solving the Lieb-Wu
equations

$\displaystyle\exp(i\kappa_{j}L)$
$\displaystyle=\prod_{\alpha=1}^{N/2}\frac{\sin\kappa_{j}-\Lambda_{\alpha}+iU/4t}{\sin\kappa_{j}-\Lambda_{\alpha}-iU/4t},$
(7)
$\displaystyle\prod_{j=1}^{N}\frac{\sin\kappa_{j}-\Lambda_{\beta}+iU/4t}{\sin\kappa_{j}-\Lambda_{\beta}-iU/4t}$
$\displaystyle=-\prod_{\alpha=1}^{N/2}\frac{\Lambda_{\alpha}-\Lambda_{\beta}+iU/2t}{\Lambda_{\alpha}-\Lambda_{\beta}-iU/2t}.$
(8)

Solving the Lieb-Wu equations via root finding can be done with great accuracy
and polynomial effort with particle number. The Bethe ansatz thus provides us
with an exact reference for the one-dimensional Hubbard chain.

Figure 1: Comparison of yrast state excitation energies for a 1D Hubbard chain
of length $L=21$ with $N=10$ fermions, obtained for interaction strengths
$U/t=0,-1,-2,-3,-4,-5$ using the Bethe ansatz (empty symbols, dashed lines)
and FCIQMC (filled symbols, solid lines). For $|U/t|\leq 5$ the agreement is
sufficient for reproducing the main features of the yrast dispersion. The
lines are merely a guide to the eye.

This allows us to compare QMC results with exact results to determine the
parameter range where we can consider FCIQMC to be reliable in the sense that
a possible initiator bias is smaller than the statistical uncertainty inherent
in the Monte Carlo approach. In general, larger values of $|U/t|$ lead to
stronger correlations in the many-body wave function, which then requires a
larger number of Slater determinants to be accurately represented. Also, the
fermionic sign problem becomes more severe, which tends to increase the
initiator bias in the calculated energy. In FIG. 1, we compare results
obtained using the Bethe ansatz and FCIQMC results for a 1D chain ($L=21$,
$W=1$). We see that for $|U/t|\leq 5$ the agreement is very good. We therefore
mainly use interaction strengths of $|U/t|\leq 5$ in this paper, which covers
the entire BEC-BCS crossover and typical values achievable in experiments
[59].

##  V Umklapp points

Figure 1 illustrates the characteristic shape of the yrast dispersion, which
is concave downward resembling inverted parabolas in the intervals
$0<P<P_{0}N/2$ and $P_{0}N/2<P<P_{0}N$, and with local minima at integer
multiples of $P_{0}N/2$. It can be understood by looking at the non-
interacting case: To increase the total momentum of the system by a single
unit of $P_{0}$, first a fermion at the Fermi surface is excited. The
resulting hole can then be filled by another fermion to increase the momentum
again, and so forth, with the energy tracing the inverted parabolic part of
the Hubbard lattice dispersion relation of Eq. (2). The first local minimum of
the yrast spectrum, called the half-umklapp point, in our case ($N=10$) at
$P/P_{0}=5$ is reached when all particles of one spin component have each been
boosted by $P_{0}$. At that point, the Fermi surfaces of both components are
shifted with respect to each other, but there are no holes in the Fermi seas
of either spin species.

The full umklapp point at $P/P_{0}=10$ is reached when both spin components,
or all particles are boosted. In the continuum limit, where full Galilean
invariance is restored, the excitation energy of the umklapp point is
determined by the boost only and is independent of interactions, since the
state is strictly a boosted ground state. In the lattice system, where
Galilean invariance is broken, a weak interaction dependence at the umklapp
points is nevertheless observed as seen in Fig. 1.

Figure 2: Yrast spectrum for a lattice of $21\times 11$ sites and interaction
values $U/t=0,-1,-2,-3,-4,-5$. The spectrum changes from having linear
segments to a parabolic shape like in the 1D case depicted in Fig. 1.

For a non-interacting 2D system, constructing the yrast dispersion from hole
excitations leads to a different shape, which in the thermodynamic limit in an
isotropic 2D system is linear. In our case, as is shown in FIG. 2 for $W=11$,
due to finite-size effects in our mesoscopic system, the spectrum for $U=0$
has linear segments but is not perfectly linear. It is remarkable that for
increasing interaction strength, the parabolic shape of the 1D spectrum with
the umklapp points at $P/P_{0}=5,10$ is restored.

Figure 3: Opposite-spin pair correlation function
$g_{\uparrow\downarrow}(k_{1},k_{2})=\langle
c_{k_{1},\uparrow}^{\dagger}c_{k_{2},\downarrow}^{\dagger}c_{k_{2},\downarrow}c_{k_{1},\uparrow}\rangle$
for the ground state (left) and the half-umklapp point (right) at $U/t=-1$
(top) and $U/t=-5$ (bottom). The crossover from a weakly interacting Fermi gas
to a BEC of bound pairs is clearly visible, as the half-umklapp point changes
its characteristics from one spin component shifted in momentum space with
respect to the other to a translation of the whole system. The half-umklapp
point of the Fermi system becomes the first full umklapp point of a fully
paired superfluid with 5 bosonic pairs.

We can take a look at how the pairing in the system changes the
characteristics of the umklapp points by calculating two-body correlation
functions. It is possible to obtain the reduced two-body density matrix

$\Gamma_{s_{1},s_{2},s_{3},s_{4}}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{3},\vec{k}_{4})=\langle\Psi_{1}|c_{\vec{k}_{1},s_{1}}^{\dagger}c_{\vec{k}_{2},s_{2}}^{\dagger}c_{\vec{k}_{3},s_{3}}c_{\vec{k}_{4},s_{4}}|\Psi_{2}\rangle$
(9)

from FCIQMC by simultaneously running two statistically independent QMC
simulations with solutions $\Psi_{1},\Psi_{2}$, to avoid biases [75]. This is
valuable even in the 1D case as obtaining the same quantity from the Bethe-
ansatz solution is not feasible. To illustrate the BEC-BCS crossover, we show
the opposite-spin pair correlation function
$g_{\uparrow\downarrow}(\vec{k}_{1},\vec{k}_{2})=\Gamma_{\uparrow\downarrow\downarrow\uparrow}(\vec{k}_{1},\vec{k}_{2},\vec{k}_{2},\vec{k}_{1})$
for the ground state and half-umklapp point in the 1D case in FIG. 3.

Strong pair correlations with $k_{1}+k_{2}=0$ clearly emerge as interaction
strength is increased from $U/t=-1$ to $U/t=-5$. The $P/P_{0}=5$ half-umklapp
point at weak interactions exhibits mostly the physics of a non-interacting
system, with two Fermi seas displaced with respect to each other. This is in
stark contrast to the situation at $U/t=-5$, where the correlation function is
the same as for $P=0$ but shifted in momentum space. The mere translation of
the pair correlation function is consistent with interpreting the system as
superfluid of 5 bosonic pairs, where total momentum $P/P_{0}=5$ correponds to
a full umklapp (i.e. Galilean boost of the ground state) in contrast to the
weakly-interacting Fermi gas, which only reaches a half umklapp point at this
momentum.

Figure 4: Excitation energy for the half umklapp point $P_{h}=5P_{0}$ (purple
data set at lower energies) and for the full umklapp point $P_{u}=10P_{0}$
(green data set at higher energies) for different values of interaction
strength $U$ and width $W$. Solid lines are calculated for $W=1$ using the
Bethe ansatz and symbols denote FCIQMC results, for the 1D system in panel (a)
and for 2D systems in panel (b). The horizontal dashed lines show the behavior
expected in the continuum 1D Fermi gas (YG model) as explained in the main
text while the dotted lines show the expected asymptotic behavior for $-U/t\gg
1$. While the excitation energy at the full umklapp point appears to be
largely independent of the transverse lattice size indicating a full
translation of both Fermi seas in momentum space, the half umklapp point shows
additional features of a mesoscopic system for different lattice sizes.
Interestingly, the values for $W=11$ approach the other points from below. For
this particular value of $W$, a rearrangement rather than a simple shift of
the fermions in momentum space takes place.

The yrast excitation energy at the half umklapp point $P_{h}=P_{0}N/2$ and at
the full umklapp point $P_{u}=P_{0}N$ are shown in Fig. 4. The lattice results
can be compared to the expected excitation energies in an equivalent free-
space system. The behavior of the umklapp points of an attractive 1D Fermi gas
in free space has been studied using the Yang-Gaudin (YG) model [47]. There,
the energy of the full umklapp point is independent of the interaction
strength, as it represents simply a translation of the entire Fermi sea in
momentum space. In the YG model this energy is given by
$E-E_{0}=P_{u}^{2}/2Nm$, where $N$ is the total particle number and
$P_{u}=P_{0}N$ the full umklapp momentum. In units of the Hubbard model
parameters, this excitation energy is $E-E_{0}=40\pi^{2}t/L^{2}$ and is
depicted as the upper dashed horizontal line in Fig. 4 (a). The half-umklapp
point however drops by factor of $2$ from $E-E_{0}=P_{h}^{2}/2mN_{\uparrow}$
to $E-E_{0}=P_{h}^{2}/2mN$, as it changes from being the half-umklapp point of
a system of $N$ fermions to being the full umklapp point of a gas of $N/2$
bosons. This is shown as the two lower dashed horizontal lines in Fig. 4 (a),
with energies $E-E_{0}=20\pi^{2}t/L^{2}$ and $E-E_{0}=10\pi^{2}t/L^{2}$.

For the non-interacting case we observe expected behavior with energy values
slightly lower than the YG model. This is because the YG model uses a
quadratic dispersion written in parameters of the Hubbard model as

$\epsilon_{\vec{k}}^{\text{YG}}=2t\alpha^{2}{\vec{k}}^{2},$ (10)

and $\epsilon_{\vec{k}}\leq\epsilon_{\vec{k}}^{\text{YG}}$.

However, for finite interactions the energy values we calculate for the
Hubbard model do not behave like for the YG model. In our system, lattice
effects dominate once the interaction is strong enough. It is known that for
$U/t\rightarrow-\infty$, the asymptotic effective Hamiltonian of the Hubbard
model corresponds to a bosonic system with a one-boson-per-site hard-core
condition and repulsive next-nearest neighbor interactions [76]. The effective
Hamiltonian also has a global pre-factor of $U^{2}/t$, meaning that the entire
spectrum has the same scaling in the asymptotic regime. For the (half-)umklapp
points, the asymptotes $80P_{0}^{2}/m\times t^{2}/U$ and $20P_{0}^{2}/m\times
t^{2}/U$ are presented in Fig. 4 (a) as the green and purple dashed lines,
respectively. We see that finite-size effects are reduced as the umklapp
energies approach these asymptotic lines.

##  VI Maxima of the yrast spectrum

The point $P=P_{h}/2$ near the first local maximum of the yrast spectrum (such
as at $P/P_{0}=2,3$ in Fig. 2) is of particular interest as a point where in
the 1D homogenous case dark solitons appear [27, 47, 77] that are stationary
with respect to background and with phase step $\pi$ across the soliton. Dark
solitons in a Fermi superfluid are characterized by a localized density
depression and a phase jump in the superfluid order parameter around this
depression. In a system with periodic boundary conditions, this phase jump
must be compensated by a phase gradient along the system, which corresponds to
a constant counterflow velocity $v_{\text{cf}}$. In addition, the soliton can
be associated with an inertial mass $m_{I}$, related to the curvature of the
yrast dispersion. We extract these parameters from our calculated dispersions
by fitting the quadratic function

$E(P)=E(0)+v_{\text{cf}}(P-P_{h}/2)+\frac{1}{2m_{I}}(P-P_{h}/2)^{2}$ (11)

around the first local maximum at momenta $P/P_{0}=1,2,3,4$, where
$P_{h}=NP_{0}/2=5P_{0}$.

Figure 5: Inertial mass of the first yrast maximum vs. interaction strength.
Unlike in free space, in the lattice, the inertial mass asymptotically
approaches $-\infty$ linearly with the interaction strength $U/t$. In the 2D
regime, where a parabolic yrast dispersion reappears for strong interactions,
the inertial mass is larger by a factor of 2.

The results for the inertial mass are shown in Fig. 5, where we only show data
points for the parameters where the yrast dispersion closely resembles a
parabolic shape. This would correspond to a regime where a superfluid is
present and the particles are strongly paired. This is mostly the case for
$W\leq 7$, where the system is still effectively one-dimensional and
transverse momentum states are sparsely populated. For the larger and more 2D
systems with $W=7,9$, we find that for $|U/t|=4,5$, a parabolic yrast spectrum
reappears (see also Fig. 2). For the effective mass, there is an increase in
magnitude by a factor of 2 as we increase the system size to $W=11$. This
indicates a further change to the properties of the system, likely a
transition from a soliton state to a vortex pair. This scenario is closely
related to the snaking instability of a planar dark soliton in a two-
dimensional superfluid, where the soliton decays into pairs of oppositely
charged vortices as the system becomes wide enough [68, 41, 70]. However, we
cannot directly show a potential density depression caused by the soliton. The
reason is that, unlike in mean-field theory, our technique provides the many-
body wave function of a (translationally invariant) eigenstate of total
momentum and thus a superposition of solitons at all possible positions. The
real-space single-particle density we can calculate is flat. To map out the
dark soliton as described in [77], we would need access to higher-order
density matrices beyond the two-body density matrix. Therefore, it remains to
be seen if the increase in inertial mass really corresponds to a soliton-
vortex pair transition.

Figure 6: Pair densities for 2D lattices with $21\times 11$ sites. The weakly-
interacting case ($U/t=-1$, top row) show a small peak in the pair density at
zero momentum both for the ground state (top left) and the yrast state at
$P/P_{0}=2$ (top right). For strong interactions ($U/t=-5$, bottom row), the
peak at zero momentum dominates at $P=0$ (bottom left), while at $P/P_{0}=2$
(bottom right), a second peak appears, with the system showing condensation
into both zero momentum and momentum $\hbar k_{x}=2P_{0}$.

Instead, we investigate more closely the pair condensation, for which a
relevant quantity is the pair Green’s function

$G_{p}(l)=\langle\psi_{j+l,\uparrow}^{\dagger}\psi_{j+l,\downarrow}^{\dagger}\psi_{j,\downarrow}\psi_{j,\uparrow}\rangle,$
(12)

where $\psi$ and $\psi^{\dagger}$ denote creation and annihilation operators,
respectively, in position space. The Fourier transform of this Green’s
function is the momentum-space pair density. For a system with a homogeneous
density, it can be directly obtained from the momentum representation of the
two-body density matrix

$n_{p}(\vec{k})=\sum_{\vec{k}_{1},\vec{k}_{2}}\Gamma_{\uparrow\downarrow\downarrow\uparrow}(\vec{k}_{1},\vec{k}-\vec{k}_{1},\vec{k}-\vec{k}_{2},\vec{k}_{2}).$
(13)

This quantity indicates whether Bose-Einstein condensation of pairs occurs.
While we find evidence of Bose-Einstein condensation of pairs by a peak in the
pair density that grows with increasing interaction strength for the ground
state and umklapp point, the situation is more complex for general yrast
states. In Fig. 6 we show the pair density for a 2D system with $21\times 11$
sites for interaction parameters $U/t=-1$ and $U/t=-5$, for the ground state
and the $P/P_{0}=2$ yrast state. For weak interactions, the structure of the
pair density is determined mostly by the structure of the Fermi sea, or the
noninteracting yrast state. For stronger interactions, the ground state
exhibits one sharp peak at zero momentum indicating strong pairing
correlations, as expected for crossover to a BEC of pairs. However, for
$P/P_{0}=2$, there are actually two peaks, for longitudinal momenta $0$ and
$2P_{0}$. Similar features appear in the smaller 2D and 1D systems.

We now show that this double-peak feature signifies the presence of fragmented
condensation. Fragmentation occurs when during the transition to a Bose-
Einstein condensate, more than one state becomes macroscopically occupied
[78]. In Fig. 6, the momentum states $(0,0)$ and $(2P_{0},0)$ dominate the
pair density. For the 1D system, where obtaining the reduced density matrices
is easier, we plot the pair densities of several momenta for the yrast state
with $P=3P_{0}$ in Fig. 7 as a function of interaction strength. We see that
in addition to $P=0$, the pair density at $P=P_{0}$ strongly increases as
well. Similarly, for other yrast states we see the same phenomenon, a strong
signature of fragmented condensation. It is worth noting that the pair density
with two peaks obtained here is similar to the case of an imbalanced Fermi
gas, where Fermi surfaces of different size lead to FFLO pairing with non-zero
total momentum and the signature is a two-peaked pair-density. This has been
studied for 1D and 2D Hubbard models [57, 44, 58, 79, 80]. Yrast states in our
balanced system start with holes in one of the Fermi seas for weak
interaction, which also leads to pairing with non-zero total momentum.

Figure 7: Pair densities for the 1D ground state (top left), the yrast states
with $P/P_{0}=2$ (top right), $P/P_{0}=3$ (bottom left) and the half-umklapp
point (bottom right). The ground state shows a rapid growth of the pair
density at zero momentum, as expected for Bose condensation of pairs, while
the half-umklapp point is identical to the ground state but shifted by one
lattice momentum unit. For the yrast states however, we observe near equal
growth of both pair momenta $P=0$ and $P=P_{0}$, meaning that fragmented
condensation into zero and non-zero momentum states is taking place.

##  VII Conclusions

In this paper, we have used the FCIQMC method to study the crossover from
weakly interacting fermions to a condensate of bosonic pairs for yrast states
in mesoscopic Fermi systems. With this method we can treat larger systems than
are accessible to the previously used deterministic CI or exact
diagonalization methods [48, 49] and can probe the full transition from a 1D
chain not just to quasi-1D rings, but also to full 2D systems. We obtain
energy spectra and reduced two-body density matrices for yrast states in the
attractive Hubbard model in these geometries.

Comparisons with exact Bethe ansatz results in 1D show very good agreement and
demonstrate that FCIQMC can accurately provide yrast states for mesoscopic
systems with 10 fermions and $21\times W$ sites where $W$ ranges from 1 to 11.
We find that the shape of the yrast spectrum changes from the typical inverted
parabolas in 1D to a quasi-linear spectrum in 2D. However, as interaction
strength is increased, the parabolic dispersion including the half umklapp
point is restored. While the quasi-linear spectrum is a consequence of the 2D
geometry and the Pauli exclusion principle, with stronger interactions the
exact shape of the non-interacting Fermi sea plays less of a role until we see
the expected universal concave downward dispersion of a spinless superfluid.
This indicates a transition to a fully paired Fermi superfluid.

We further find that mesoscopic effects can cause significant deviations for
certain geometries from the general behavior of the umklapp energies, which
otherwise does not differ much between 1D and 2D systems. Specifically we find
that the half-umklapp excitation energy for a lattice of $21\times 11$ points
increases with interaction strength at intermediate values of $U/t$ contrary
to the general trend displayed by all other systems under study. This
originates in a rearrangement of the fermions in momentum space, where for
this particular geometry, the half-umklapp point is a different configuration
than the ground state with one spin component shifted in momentum space. This
can be of importance for experiments on mesoscopic Fermi systems.

By calculating the pair correlation function, we can follow the pairing
process by means of which the fermionic half-umklapp point becomes the first
full umklapp point of the Bose condensate of pairs.

In the fully paired regime, we calculated the inertial mass of the dark-
soliton-like local maximum of the yrast dispersion. We found a sudden increase
by a factor of 2 between the narrow stripe geometry and our largest system in
a lattice of $21\times 11$ sites. We interpret this as a possible change in
the system geometry where a soliton is no longer stable and the yrast state is
instead provided by pair of oppositely charged vortices.

The most striking feature of yrast states in the Hubbard model is revealed to
be fragmented condensation, which occurs around the maxima of the yrast
dispersion, away from the umklapp points. We find that here, multi-peaked pair
densities appear, where in addition to pairing with zero total momentum, the
amplitude of other total momentum pairs becomes large for strong interactions.
These pair densities of yrast states share some similarities with FFLO states,
which are characterized by a double-peaked momentum pair density [44]. In the
FFLO case the origin of this is a mismatch of Fermi surfaces, which have
different size due to the spin imbalance. In our case of the yrast states, the
Fermi surfaces are shifted in momentum space. In addition, both yrast and FFLO
states are related to solitons in the real space density [44].

###### Acknowledgements.

We thank M. Zwierlein, S. Shamailov, P. Jeszenszki and W. Dobrautz for useful
discussions. UE thanks the Max-Planck-Institute for Solid State Research for
hospitality. This work was supported by the Marsden Fund of New Zealand
(Contract No. MAU1604), from government funding managed by the Royal Society
Te Apārangi. We also acknowledge support by the NeSI high-performance
computing facilities through a Merit allocation and a consultancy project.

## References

  * Eagles [1969] DM Eagles, “Possible Pairing without Superconductivity at Low Carrier Concentrations in Bulk and Thin-Film Superconducting Semiconductors,” Phys. Rev. 186, 456–463 (1969).
  * Leggett [1980] Anthony J Leggett, “Diatomic Molecules and Cooper pairs,” in  _Mod. Trends Theory Condens. Matter_, Lecture Notes in Physics, edited by Andrzej Pȩkalski and JerzyA Przystawa (Springer Berlin Heidelberg, 1980) pp. 13–27.
  * Zwerger [2012] Wilhelm Zwerger, ed.,  _The BCS –- BEC Crossover and the Unitary Fermi Gas_, Lecture Notes in Physics, Vol. 836 (Springer Berlin / Heidelberg,, Berlin, Heidelberg, 2012) p. 532, arXiv:1008.3933 .
  * Levinsen and Parish [2015] Jesper Levinsen and Meera M. Parish, “Strongly interacting two-dimensional Fermi gases,” in  _Annu. Rev. Cold Atoms Mol._ (WORLD SCIENTIFIC, 2015) pp. 1–75.
  * Giorgini and Stringari [2008] Stefano Giorgini and Sandro Stringari, “Theory of ultracold atomic Fermi gases,” Rev. Mod. Phys. 80, 1215–1274 (2008), arXiv:0706.3360 .
  * Parish [2014] Meera M. Parish, “The BCS–BEC Crossover,” in  _Quantum Gas Exp._, edited by Päivi Törmä and Klaus Sengstock (Imperial College Press, 2014) pp. 179–197, arXiv:1402.5171 .
  * Regal  _et al._ [2004] C. A. Regal, M. Greiner, and D. S. Jin, “Observation of Resonance Condensation of Fermionic Atom Pairs,” Phys. Rev. Lett. 92, 040403 (2004), arXiv:0401554 [cond-mat] .
  * Zwierlein  _et al._ [2004] M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. Raupach, A. J. Kerman, and W. Ketterle, “Condensation of Pairs of Fermionic Atoms near a Feshbach Resonance,” Phys. Rev. Lett. 92, 120403 (2004).
  * Chin  _et al._ [2004] C. Chin, M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, J. Hecker Denschlag, and R. Grimm, “Observation of the Pairing Gap in a Strongly Interacting Fermi Gas,” Science 305, 1128–1130 (2004), arXiv:0405632 [cond-mat] .
  * Bourdel  _et al._ [2004] T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, F. Chevy, M. Teichmann, L. Tarruell, S. J. J. M. F. Kokkelmans, and C. Salomon, “Experimental Study of the BEC-BCS Crossover Region in Lithium 6,” Phys. Rev. Lett. 93, 050401 (2004), arXiv:0403091 [cond-mat] .
  * Kinast  _et al._ [2004] J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, and J. E. Thomas, “Evidence for Superfluidity in a Resonantly Interacting Fermi Gas,” Phys. Rev. Lett. 92, 150402 (2004).
  * Ketterle and Zwierlein [2008] Wolfgang Ketterle and Martin W Zwierlein, “Making, probing and understanding ultracold Fermi gases,” Riv. del Nuovo Cim. 31, 247 (2008), arXiv:0801.2500 .
  * Ku  _et al._ [2012] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W. Zwierlein, “Revealing the Superfluid Lambda Transition in the Universal Thermodynamics of a Unitary Fermi Gas,” Science 335, 563–567 (2012).
  * Carcy  _et al._ [2019] C. Carcy, S. Hoinka, M. G. Lingham, P. Dyke, C. C. N. Kuhn, H. Hu, and C. J. Vale, “Contact and Sum Rules in a Near-Uniform Fermi Gas at Unitarity,” Phys. Rev. Lett. 122, 203401 (2019), arXiv:1902.07853 .
  * Mukherjee  _et al._ [2019] Biswaroop Mukherjee, Parth B. Patel, Zhenjie Yan, Richard J. Fletcher, Julian Struck, and Martin W. Zwierlein, “Spectral Response and Contact of the Unitary Fermi Gas,” Phys. Rev. Lett. 122, 203402 (2019), arXiv:1902.08548 .
  * Bakr  _et al._ [2009] Waseem S Bakr, Jonathon I Gillen, Amy Peng, Simon Folling, and Markus Greiner, “A quantum gas microscope for detecting single atoms in a Hubbard-regime optical lattice,” Nature 462, 74–77 (2009).
  * Cheuk  _et al._ [2015] Lawrence W Cheuk, Matthew A Nichols, Melih Okan, Thomas Gersdorf, Vinay V Ramasesh, Waseem S Bakr, Thomas Lompe, and Martin W Zwierlein, “Quantum-Gas Microscope for Fermionic Atoms,” Phys. Rev. Lett. 114, 193001 (2015).
  * Haller  _et al._ [2015] Elmar Haller, James Hudson, Andrew Kelly, Dylan A. Cotta, Bruno Peaudecerf, Graham D. Bruce, and Stefan Kuhr, “Single-atom imaging of fermions in a quantum-gas microscope,” Nat. Phys. 11, 738–742 (2015).
  * Brown  _et al._ [2020] Peter T. Brown, Elmer Guardado-Sanchez, Benjamin M. Spar, Edwin W. Huang, Thomas P. Devereaux, and Waseem S. Bakr, “Angle-resolved photoemission spectroscopy of a Fermi–Hubbard system,” Nat. Phys. 16, 26–31 (2020), arXiv:1903.05678 .
  * Grünzweig  _et al._ [2010] T. Grünzweig, A. Hilliard, M. McGovern, and M. F. Andersen, “Near-deterministic preparation of a single atom in an optical microtrap,” Nat. Phys. 6, 951–954 (2010).
  * Serwane  _et al._ [2011] F Serwane, G Zürn, T Lompe, T B Ottenstein, A N Wenz, and S Jochim, “Deterministic preparation of a tunable few-fermion system.” Science 332, 336–338 (2011), arXiv:1101.2124 .
  * Wenz  _et al._ [2013] A Wenz, G Zürn, S Murmann, I Brouzos, T Lompe, and S Jochim, “From Few to Many: Observing the Formation of a Fermi Sea One Atom at a Time,” Science 342, 457–460 (2013), arXiv:arXiv:1307.3443v2 .
  * Reynolds  _et al._ [2020] L. A. Reynolds, E. Schwartz, U. Ebling, M. Weyland, J. Brand, and M. F. Andersen, “Direct Measurements of Collisional Dynamics in Cold Atom Triads,” Phys. Rev. Lett. 124, 073401 (2020), arXiv:2001.05141 .
  * Kulish  _et al._ [1976] P. P. Kulish, S. V. Manakov, and L. D. Faddeev, “Comparison of the exact quantum and quasiclassical results for a nonlinear Schrödinger equation,” Theor. Math. Phys. 28, 615–620 (1976).
  * Kanamoto  _et al._ [2008] Rina Kanamoto, Lincoln D. Carr, and Masahito Ueda, “Topological Winding and Unwinding in Metastable Bose-Einstein Condensates,” Phys. Rev. Lett. 100, 060401 (2008).
  * Kanamoto  _et al._ [2010] R. Kanamoto, L. D. Carr, and M. Ueda, “Metastable quantum phase transitions in a periodic one-dimensional Bose gas. II. Many-body theory,” Phys. Rev. A 81, 023625 (2010), arXiv:0910.2805 .
  * Jackson  _et al._ [2011] A. D. Jackson, J. Smyrnakis, M. Magiropoulos, and G. M. Kavoulakis, “Solitary waves and yrast states in bose-einstein condensed gases of atoms,” EPL (Europhysics Letters) 95, 30002 (2011).
  * Fialko  _et al._ [2012] Oleksandr Fialko, Marie-Coralie Delattre, Joachim Brand, and Andrey R. Kolovsky, “Nucleation in finite topological systems during continuous metastable quantum phase transitions,” Phys. Rev. Lett. 108, 250402 (2012).
  * Sato  _et al._ [2012] Jun Sato, Rina Kanamoto, Eriko Kaminishi, and Tetsuo Deguchi, “Exact relaxation dynamics of a localized many-body state in the 1D Bose gas,” Phys. Rev. Lett. 108, 110401 (2012), arXiv:1112.4244 .
  * Syrwid and Sacha [2015] Andrzej Syrwid and Krzysztof Sacha, “Lieb-Liniger model: Emergence of dark solitons in the course of measurements of particle positions,” Phys. Rev. A 92, 032110 (2015), arXiv:1505.06586 .
  * Shamailov and Brand [2019] Sophie S. Shamailov and Joachim Brand, “Quantum dark solitons in the one-dimensional bose gas,” Phys. Rev. A 99, 043632 (2019).
  * Tsuzuki [1971] Toshio Tsuzuki, “Nonlinear waves in the Pitaevskii-Gross equation,” J. Low Temp. Phys. 4, 441 (1971).
  * Yefsah  _et al._ [2013] Tarik Yefsah, Ariel T Sommer, Mark J H Ku, Lawrence W. Cheuk, Wenjie Ji, Waseem S Bakr, and Martin W Zwierlein, “Heavy solitons in a fermionic superfluid.” Nature 499, 426–30 (2013), arXiv:1302.4736 .
  * Ku  _et al._ [2014] Mark J. H. Ku, Wenjie Ji, Biswaroop Mukherjee, Elmer Guardado-Sanchez, Lawrence W Cheuk, Tarik Yefsah, and Martin W Zwierlein, “Motion of a Solitonic Vortex in the BEC-BCS Crossover,” Phys. Rev. Lett. 113, 065301 (2014), arXiv:1402.7052 .
  * Ku  _et al._ [2016] Mark J. H. Ku, Biswaroop Mukherjee, Tarik Yefsah, and Martin W Zwierlein, “Cascade of Solitonic Excitations in a Superfluid Fermi gas: From Planar Solitons to Vortex Rings and Lines,” Phys. Rev. Lett. 116, 045304 (2016), arXiv:1507.01047 .
  * Antezza  _et al._ [2007] Mauro Antezza, Franco Dalfovo, Lev Pitaevskii, and Sandro Stringari, “Dark solitons in a superfluid Fermi gas,” Phys. Rev. A 76, 043610 (2007), arXiv:0706.0601 .
  * Liao and Brand [2011] Renyuan Liao and Joachim Brand, “Traveling dark solitons in superfluid Fermi gases,” Phys. Rev. A 83, 041604(R) (2011).
  * Scott  _et al._ [2011] R. Scott, F. Dalfovo, L. Pitaevskii, and S. Stringari, “Dynamics of Dark Solitons in a Trapped Superfluid Fermi Gas,” Phys. Rev. Lett. 106, 185301 (2011).
  * Spuntarelli  _et al._ [2011] Andrea Spuntarelli, Lincoln D Carr, Pierbiagio Pieri, and Giancarlo C Strinati, “Gray solitons in a strongly interacting superfluid Fermi gas,” New J. Phys. 13, 035010 (2011).
  * Scott  _et al._ [2012] R. G. Scott, F. Dalfovo, L. P. Pitaevskii, S. Stringari, O. Fialko, R. Liao, and J. Brand, “The decay and collisions of dark solitons in superfluid Fermi gases,” New J. Phys. 14, 023044 (2012), arXiv:1109.6444 .
  * Cetoli  _et al._ [2013] A Cetoli, J Brand, R G Scott, F Dalfovo, and L P Pitaevskii, “Snake instability of dark solitons in fermionic superfluids,” Phys. Rev. A 88, 043639 (2013), arXiv:1307.3717 .
  * Efimkin and Galitski [2015] Dmitry K. Efimkin and Victor Galitski, “Moving solitons in a one-dimensional fermionic superfluid,” Phys. Rev. A 91, 023616 (2015), arXiv:1408.6511 .
  * Yoshida and Yip [2007] Nobukatsu Yoshida and S.-K. Yip, “Larkin-Ovchinnikov state in resonant Fermi gas,” Phys. Rev. A 75, 063601 (2007).
  * Lutchyn  _et al._ [2011] Roman M. Lutchyn, Maxim Dzero, and Victor M. Yakovenko, “Spectroscopy of the soliton lattice formation in quasi-one-dimensional fermionic superfluids with population imbalance,” Phys. Rev. A 84, 033609 (2011).
  * Fulde and Ferrell [1964] Peter Fulde and Richard A. Ferrell, “Superconductivity in a Strong Spin-Exchange Field,” Phys. Rev. 135, A550–A563 (1964).
  * Larkin and Ovchinnikov [1965] A.I. Larkin and Yu.N. Ovchinnikov, “Nonuniform state of superconductors,” Sov. Phys. JETP 20, 762 (1965).
  * Shamailov and Brand [2016] Sophie S Shamailov and Joachim Brand, “Dark-soliton-like excitations in the yang–gaudin gas of attractively interacting fermions,” New Journal of Physics 18, 075004 (2016).
  * Viefers  _et al._ [2004] S. Viefers, P. Koskinen, P. Singha Deo, and M. Manninen, “Quantum rings for beginners: energy spectra and persistent currents,” Phys. E Low-dimensional Syst. Nanostructures 21, 1–35 (2004), arXiv:0310064 [cond-mat] .
  * Manninen  _et al._ [2012] M Manninen, S Viefers, and S.M. Reimann, “Quantum rings for beginners II: Bosons versus fermions,” Phys. E Low-dimensional Syst. Nanostructures 46, 119–132 (2012).
  * Fomin [2018] Vladimir M. Fomin, ed.,  _Physics of Quantum Rings_, 2nd ed., NanoScience and Technology (Springer International Publishing, Cham, 2018) p. 585.
  * Troyer and Wiese [2005] Matthias Troyer and Uwe-Jens Wiese, “Computational Complexity and Fundamental Limitations to Fermionic Quantum Monte Carlo Simulations,” Phys. Rev. Lett. 94, 170201 (2005), arXiv:0408370 [cond-mat] .
  * Wu and Zhang [2005] Congjun Wu and Shou-Cheng Zhang, “Sufficient condition for absence of the sign problem in the fermionic quantum monte carlo algorithm,” Phys. Rev. B 71, 155115 (2005).
  * Foulkes  _et al._ [2001] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, “Quantum monte carlo simulations of solids,” Rev. Mod. Phys. 73, 33–83 (2001).
  * Carlson  _et al._ [2011] J. Carlson, Stefano Gandolfi, Kevin E. Schmidt, and Shiwei Zhang, “Auxiliary-field quantum monte carlo method for strongly paired fermions,” Phys. Rev. A 84, 061602 (2011).
  * Motta  _et al._ [2019] Mario Motta, Shiwei Zhang, and Garnet Kin-Lic Chan, “Hamiltonian symmetries in auxiliary-field quantum monte carlo calculations for electronic structure,” Phys. Rev. B 100, 045127 (2019).
  * Varney  _et al._ [2009] C. N. Varney, C.-R. Lee, Z. J. Bai, S. Chiesa, M. Jarrell, and R. T. Scalettar, “Quantum monte carlo study of the two-dimensional fermion hubbard model,” Phys. Rev. B 80, 075116 (2009).
  * Wolak  _et al._ [2010] M. J. Wolak, V. G. Rousseau, C. Miniatura, B. Grémaud, R. T. Scalettar, and G. G. Batrouni, “Finite-temperature quantum monte carlo study of the one-dimensional polarized fermi gas,” Phys. Rev. A 82, 013614 (2010).
  * Wolak  _et al._ [2012] M. J. Wolak, B. Grémaud, R. T. Scalettar, and G. G. Batrouni, “Pairing in a two-dimensional fermi gas with population imbalance,” Phys. Rev. A 86, 023630 (2012).
  * Mitra  _et al._ [2018] Debayan Mitra, Peter T. Brown, Elmer Guardado-Sanchez, Stanimir S. Kondov, Trithep Devakul, David A. Huse, Peter Schauß, and Waseem S. Bakr, “Quantum gas microscopy of an attractive Fermi–Hubbard system,” Nat. Phys. 14, 173–177 (2018), arXiv:1705.02039 .
  * Booth  _et al._ [2009] George H. Booth, Alex J. W. Thom, and Ali Alavi, “Fermion monte carlo without fixed nodes: A game of life, death, and annihilation in slater determinant space,” The Journal of Chemical Physics 131, 054106 (2009), https://aip.scitation.org/doi/pdf/10.1063/1.3193710 .
  * Cleland  _et al._ [2010] Deidre Cleland, George H. Booth, and Ali Alavi, “Communications: Survival of the fittest: Accelerating convergence in full configuration-interaction quantum monte carlo,” The Journal of Chemical Physics 132, 041103 (2010), https://doi.org/10.1063/1.3302277 .
  * Booth  _et al._ [2013] George H Booth, Andreas Grüneis, Georg Kresse, and Ali Alavi, “Towards an exact description of electronic wavefunctions in real solids.” Nature 493, 365–70 (2013).
  * Cleland  _et al._ [2012] Deidre Cleland, George H Booth, Catherine Overy, and Ali Alavi, “Taming the first-row diatomics: A full configuration interaction quantum Monte Carlo study,” J. Chem. Theory Comput. 8, 4138–4152 (2012).
  * Jeszenszki  _et al._ [2020] Péter Jeszenszki, Ulrich Ebling, Hongjun Luo, Ali Alavi, and Joachim Brand, “Eliminating the wave-function singularity for ultracold atoms by a similarity transformation,” Phys. Rev. Res. 2, 043270 (2020), arXiv:2002.05987 .
  * Yang  _et al._ [2020] Mingrui Yang, Elke Pahl, and Joachim Brand, “Improved walker population control for full configuration interaction quantum Monte Carlo,” J. Chem. Phys. 153, 174103 (2020), arXiv:2008.01927 .
  * Spencer  _et al._ [2012] J. S. Spencer, N. S. Blunt, and W. M.C. Foulkes, “The sign problem and population dynamics in the full configuration interaction quantum monte carlo method,” The Journal of Chemical Physics 136, 054110 (2012), https://doi.org/10.1063/1.3681396 .
  * Brand and Reinhardt [2001] Joachim Brand and William P Reinhardt, “Generating ring currents, solitons and svortices by stirring a Bose-Einstein condensate in a toroidal trap,” J. Phys. B At. Mol. Opt. Phys. 34, L113–L119 (2001).
  * Brand and Reinhardt [2002] Joachim Brand and William P. Reinhardt, “Solitonic vortices and the fundamental modes of the ”snake instability”: Possibility of observation in the gaseous Bose-Einstein condensate,” Phys. Rev. A 65, 043612 (2002).
  * Muñoz Mateo and Brand [2014] A. Muñoz Mateo and J. Brand, “Chladni Solitons and the Onset of the Snaking Instability for Dark Solitons in Confined Superfluids,” Phys. Rev. Lett. 113, 255302 (2014), arXiv:1408.0947 .
  * Van Alphen  _et al._ [2019] W. Van Alphen, H. Takeuchi, and J. Tempere, “Crossover between snake instability and josephson instability of dark solitons in superfluid fermi gases,” Phys. Rev. A 100, 023628 (2019).
  * Muñoz Mateo and Brand [2015] A. Muñoz Mateo and J. Brand, “Stability and dispersion relations of three-dimensional solitary waves in trapped Bose-Einstein condensates,” New J. Phys. 17, 125013 (2015), arXiv:1510.01465 .
  * Toikka and Brand [2017] L. A. Toikka and J. Brand, “Asymptotically solvable model for a solitonic vortex in a compressible superfluid,” New J. Phys. 19, 023029 (2017), arXiv:1608.08701 .
  * Dobrautz  _et al._ [2019] Werner Dobrautz, Hongjun Luo, and Ali Alavi, “Compact numerical solutions to the two-dimensional repulsive hubbard model obtained via nonunitary similarity transformations,” Phys. Rev. B 99, 075119 (2019).
  * Lieb and Wu [2003] Elliott H. Lieb and F.Y. Wu, “The one-dimensional hubbard model: a reminiscence,” Physica A: Statistical Mechanics and its Applications 321, 1 – 27 (2003), statphys-Taiwan-2002: Lattice Models and Complex Systems.
  * Overy  _et al._ [2014] Catherine Overy, George H. Booth, N. S. Blunt, James J. Shepherd, Deidre Cleland, and Ali Alavi, “Unbiased reduced density matrices and electronic properties from full configuration interaction quantum Monte Carlo,” J. Chem. Phys. 141, 244117 (2014).
  * Teubel  _et al._ [1990] A Teubel, E Kolley, and W Kolley, “On a unified canonical transformation of the large-negative-(positive-) u hubbard model,” Journal of Physics A: Mathematical and General 23, L837–L840 (1990).
  * Syrwid  _et al._ [2018] Andrzej Syrwid, Dominique Delande, and Krzysztof Sacha, “Emergence of dark soliton signatures in a one-dimensional unpolarized attractive fermi gas on a ring,” Phys. Rev. A 98, 023616 (2018).
  * Mueller  _et al._ [2006] Erich J. Mueller, Tin-Lun Ho, Masahito Ueda, and Gordon Baym, “Fragmentation of bose-einstein condensates,” Phys. Rev. A 74, 033612 (2006).
  * Cheng  _et al._ [2018a] Song Cheng, Yuzhu Jiang, Yi-Cong Yu, Murray T. Batchelor, and Xi-Wen Guan, “Asymptotic correlation functions and fflo signature for the one-dimensional attractive hubbard model,” Nuclear Physics B 929, 353 – 376 (2018a).
  * Cheng  _et al._ [2018b] Song Cheng, Yi-Cong Yu, M. T. Batchelor, and Xi-Wen Guan, “Fulde-ferrell-larkin-ovchinnikov correlation and free fluids in the one-dimensional attractive hubbard model,” Phys. Rev. B 97, 121111 (2018b).


#  LLOL: Low-Latency Odometry for Spinning Lidars

Chao Qu, Shreyas S. Shivakumar, Wenxin Liu and Camillo J. Taylor We gratefully
acknowledge the support of Distributed and Collaborative Intelligent Systems
and Technology Collaborative Research Alliance (DCIST) and C-BRIC, a
Semiconductor Research Corporation Joint University Microelectronics Program,
program cosponsored by DARPA.C. Qu, S. S. Shivakumar, W. Liu and C. J. Taylor
are with the GRASP Laboratory, School of Engineering and Applied Sciences,
University of Pennsylvania {quchao, sshreyas, wenxinl, cjtaylor}
@seas.upenn.edu

###### Abstract

In this paper, we present a low-latency odometry system designed for spinning
lidars. Many existing lidar odometry methods wait for an entire sweep from the
lidar before processing the data. This introduces a large delay between the
first laser firing and its pose estimate. To reduce this latency, we treat the
spinning lidar as a streaming sensor and process packets as they arrive. This
effectively distributes expensive operations across time, resulting in a very
fast and lightweight system with a much higher throughput and lower latency.
Our open source implementation is available at
https://github.com/versatran01/llol.

##  I Introduction

Light Detection and Ranging (lidar) sensors are one of the primary sensing
modalities in autonomous driving and robotics applications. They work by
sweeping one or more laser beams through the scene and calculate range by
measuring time of flight. This is apparent in mechanical lidars that
continually spin a comb of lasers. However, solid-state lidars designed to
minimize moving parts also involve redirecting a laser beam over time to cover
a wider field of view (FOV).

Many existing lidar odometry systems wait for the lidar to complete a full
revolution and then process the resulting sweep to update ego-motion along
with any associated scene representations. Because of this, systems that
combine Inertial Measurement Units (IMUs), cameras and lidars often accept
that lidar is the slowest sensor in the system, and design algorithms based on
that assumption [1, 2, 3]. At first glance this is a valid assumption since a
spinning lidar usually rotates at $10\text{\,}\mathrm{H}\mathrm{z}$, which is
a much lower rate compared to those of a typical IMU
($>$$100\text{\,}\mathrm{H}\mathrm{z}$) or camera
($>$$20\text{\,}\mathrm{H}\mathrm{z}$). However, this data capture delay is an
abstraction of the underlying mechanism of the lidar, and is adopted only for
ease of processing. It ignores the streaming nature of the sensor and the fact
that the raw lidar packets usually arrive at a much higher rate
($>$$600\text{\,}\mathrm{H}\mathrm{z}$). Recent works have also demonstrated
that directly operating on lidar packets enables low-latency object detection,
which is critical for autonomous driving [4, 5].

Figure 1: Difference in latency when using full sweep (top) vs partial sweeps
(bottom) in an odometry system. Gray blocks denote lidar acquisition time, and
yellow blocks denote odometry runtime. We propose to process partial sweeps
which allows the system to produce odometry at a higher rate with lower
latency.

In this paper we show how to exploit this continuous stream of data to
significantly decrease latency and increase throughput for lidar odometry.
Fig. 1 shows the difference in latency of a standard lidar odometry and the
one proposed in this paper. Consider the latency between the time of the laser
firing at the beginning of the sweep and the time that the sweep trajectory is
estimated, we note that it has two components: the integration time of the
lidar and the runtime of the odometry. By this definition, a system that
operates on a full sweep will have roughly twice the latency compared to one
that operates on half a sweep.

Our proposed approach has two advantages. Firstly it can process partial
sweeps by employing a circular buffer of range measurements which drastically
increases the rate at which we can produce new pose measurements, for example
our approach can produce pose estimates at a rate of 80Hz from a lidar sensor
with a sweep rate of 10Hz. Secondly our approach can process the resulting
range information extremely efficiently on modern processor architectures
without the need for GPU acceleration which minimizes the latency of those
pose estimates. Both of these properties, increased throughput and decreased
latency, are particularly important in contexts where the pose estimates are
being used for closed loop control of dynamic platforms like quadrotors.

In the rest of this paper, we describe key elements of our system in Section
III and provide implementation details in Section IV. We then evaluate its
accuracy and runtime by comparing to a state-of-the-art lidar odometry system
[6] in Section V. We establish that by using partial sweeps we can increase
the throughput substantially while lowering the effective latency. We show
that by adopting a data-oriented design with array-based data structures, our
system exhibits almost linear scalability in terms of multi-threading
performance. Unfortunately, many popular datasets only provide lidar data pre-
packaged as point clouds from full sweeps. Therefore, we also collected our
own dataset to further validate our method. Our system, named LLOL, is open
sourced and available to the community.

##  II Related Work

In this section, we review related work on lidar-(inertial) odometry with a
particular focus on systems that have been demonstrated to run in real time on
embedded platforms.

###  II-A Lidar Odometry

In a seminal paper Zhang et al. proposed LOAM, a feature based lidar odometry
designed to work on embedded platforms, with loose IMU integration and a
discrete trajectory representation [7]. Shan et al. improved LOAM for ground
robots by leveraging an initial segmentation procedure to detect ground planes
[8]. LOAM was further extended by Lin et al. to work with limited FOV LiDARs,
supporting irregular sampling of feature points and enhanced outlier
rejection. Liu et al. proposed a method of local bundle adjustment and a
temporal sliding window that allows for real-time performance while also
improving odometry, and use LOAM as a front end in their validation [9].
Elastic-Fusion, SuMa and MARS use a surfel based map to improve performance,
and allow representation of larger scale environments [10, 11, 12, 13].
Palieri et al. proposed LOCUS, a multi-sensor lidar-centric system for
odometry and mapping that is robust to intermittent sensor degradation, and
deployed it on both legged and wheeled robot platforms. Yokozuka et al.
proposed LITAMIN [14], which is a light-weight lidar odometry. They adopt a
voxelization step to reduce the number of points used in ICP and employ a
modified GICP [15] cost that avoids expensive eigen-decomposition computations
commonly used with the point-to-plane metric. LITAMIN2 further improves real-
time performance by matching two distributions via a symmetric KL-divergence
cost [16].

###  II-B Lidar-inertial Odometry

Neuhaus et al. propose MC2SLAM, which is a tightly-coupled odometry system
that is designed with real-time performance in mind. They report a runtime of
approximately $60\text{\,}\mathrm{m}\mathrm{s}$, excluding an intermittent
pose graph optimization step that requires an additional
$108\text{\,}\mathrm{m}\mathrm{s}$ every five frames. In [17], the authors
propose LiLi-OM to address the irregular scan pattern of the Livox lidar
similar to [18]. They perform tightly coupled lidar-inertial fusion with a
keyframe-based sliding window optimization. They report a runtime performance
of roughly $85\text{\,}\mathrm{m}\mathrm{s}$ per sweep on Velodyne lidar data.
LIO-SAM [19] improves upon [20] and uses a window of local sweeps instead of
matching individual lidar sweeps to a global map. They report an improvement
in efficiency and accuracy, and a runtime performance of roughly
$50\text{\,}\mathrm{m}\mathrm{s}$-$100\text{\,}\mathrm{m}\mathrm{s}$ per sweep
on a laptop CPU. Xu et al. proposed Fast-LIO and Fast-LIO2, which we use in
our quantitative comparisons given their focus on real-time performance and
support for similar hardware as our method [21, 6]. Fast-LIO2 achieves real-
performance on embedded CPUs via the use of an incremental KD Tree and
tightly-couples IMU data. Perhaps the most similar work to ours is UPSLAM
[22], which uses a union of depth panoramas [23] as its map representation. It
is also a full-fledged SLAM system that includes loop closure. However, UPSLAM
requires a GPU since it utilizes all points from a lidar sweep.

##  III Method

A spinning lidar sensor produces data in packets which can be decoded by a
driver to produce a final point cloud. It can also be converted into a
spherical coordinate system resulting in a range image. In this work, we
define a lidar sweep to be either a point cloud or a range image generated
from a complete revolution of the sensor.

Our system consumes lidar packets directly and accumulates a partial sweep as
a range image internally (Fig. 2(a)). Compared with structured point clouds,
range images are more compact and memory efficient. Once the partial sweep
reaches a desired integration time, the odometry system processes and produces
a pose estimate. For ease of illustration, we assume a full sweep as input in
this section. We address how to adapt our system to partial sweeps in section
IV.

###  III-A Feature Detection

(a) Input sweep

(b) Cell smoothness

(c) Cell variance

(d) Match candidates

(e) Matched cells

Figure 2: Various stages in the feature detection process. Brighter color
denotes higher value in both colormaps.

Figure 3: Depth panoramas generated by our system for both indoor (top) and
outdoor (bottom) scenes using the same resolution ($256\times 1024$) with a
90° vertical FOV.

The input to our system is a range image from which we can recover the 3D
position of each pixel. There is usually too much information to process, as a
typical lidar sweep consists of $64\times 1024=65536$ data points. To reduce
the amount of information, we divide the sweep into a grid of 2D cells and
compute smoothness (2(b)) and variance (2(c)) for each cell. Similar to [8],
smoothness $s$ is defined as

$\displaystyle s=\frac{1}{|C|}\left|\sum_{n\in C}\frac{r_{n}}{r_{m}}-1\right|$
(1)

where $r_{n}$ is the range of pixel $n$, $r_{m}$ is the range of the midpoint
and $|C|$ is the number of valid pixels in cell $C$.

We then filter the entire grid by a maximum smoothness and variance threshold
to identify planar regions in the sweep. An optional non-minimum suppression
can be applied to further reduce the number of cells while maintaining a
uniform distribution across the grid. We represent each remaining cell (2(d))
with a multivariate Gaussian with mean $\mu^{\mathcal{S}}\in\mathbb{R}^{3}$
and covariance $\Sigma^{\mathcal{S}}\in\mathbb{R}^{3\times 3}$ which captures
the statistics of the 3D positions of the points in that cell. These
quantities will be used for data association and optimization. This feature
detection stage is similar to the voxelization step in LITAMIN2 [16] but is
carried out in the 2D image space.

###  III-B Local Map Representation

Following UPSLAM [22], we employ a depth panorama (pano) [23] as our local
map. A pano is similar to a sweep, but with higher resolution and represents
the 3D structure at one point along the sensors trajectory. It is built by
fusing multiple localized and motion-corrected sweeps. Compared with 3D map
representations like point clouds or voxel/surfel grids, a depth pano has a
fixed storage complexity. The advantages include low memory footprint, linear
access patterns and adaptive spatial resolution. As shown in Fig. 3, the same
pano size can be used to model both indoor and outdoor scenes. Whereas with 3D
representations, one has to choose a suitable spatial resolution to voxelize
the workspace, which results in high variance in runtime efficiency.

Due to its fixed viewpoint and spherical projection model, our pano map is
susceptible to occlusions and limited horizontal FOV. In practice, we found
this to be less of a problem since most of our robots move in open space with
planar motions. Note that a graph of pano s can be used for pose graph
optimization to reduce drift [22], or fed into other loop closure systems [24,
25].

###  III-C Local Trajectory Representation

At any given time, we maintain a local trajectory of the current sweep, which
has a time span of $[t_{0},t_{1})$. It consists of a set of discrete states,
each for a single column of cells in the grid. Given the laser firing time
$\Delta t$, we can recover the column time by $t_{c}=t_{0}+c\cdot\Delta t$.
Each state $X^{\mathcal{P}}_{\mathcal{S}_{t}}$ consists of the position,
velocity and orientation of the sensor at time $t\in[t_{0},t_{1})$ with
respect to the current pano frame

$\displaystyle
X^{\mathcal{P}}_{\mathcal{S}_{t}}=\begin{bmatrix}p^{\mathcal{P}}_{\mathcal{S}_{t}}&v^{\mathcal{P}}_{\mathcal{S}_{t}}&R^{\mathcal{P}}_{\mathcal{S}_{t}}\end{bmatrix}$
(2)

Given a new sweep, we start from the last state of the previous sweep and
propagate forward in time using the available IMU readings to obtain an
initial guess. This trajectory estimate is then used to resolve the locations
of the features in the sweep.

###  III-D Projective ICP

Figure 4: Data association between sweep and pano. Red and green ellipsoids
depict the normal distributions of grid cells and their matching windows in
the pano.

Given a set of match candidates (planar cells from III-A), we can estimate the
pose of the current sweep by aligning it to the pano via an Iterative closest
point (ICP) procedure.

####  III-D1 Data Association

We find correspondences by projective data association [26], where points in
the sweep frame $\mathcal{S}$ are transformed into the pano frame
$\mathcal{P}$

$\displaystyle
u^{\mathcal{P}}=\pi(T^{\mathcal{P}}_{\mathcal{S}}\mu^{\mathcal{S}})=\pi(R^{\mathcal{P}}_{\mathcal{S}}\mu^{\mathcal{S}}+p^{\mathcal{P}}_{\mathcal{S}}),$
(3)

Here $\mu^{\mathcal{S}}$ is the mean position of the cell from the sweep and
$u^{\mathcal{P}}$ is its target pixel coordinate in the pano. $\pi(\cdot)$ is
the spherical projection function, and
$T^{\mathcal{P}}_{\mathcal{S}}\in\mathbb{SE}(3)$ is the rigid transformation
from the sweep to the pano frame. We only accept matches that pass the
following checks: 1) the target pixel must lie within the pano, 2) the depths
of the target pixel and sweep point must be close enough, and 3) more than
half of the pixels within a window around the target must be valid. These
criteria will eliminate most of the outliers due to occlusion. The remaining
ones can be handled by either a robust loss function or $\chi$-squared test.
We summarize the target pixel with a multivariate Gaussian with mean
$\mu^{\mathcal{P}}\in\mathbb{R}^{3}$ and covariance
$\Sigma^{\mathcal{P}}\in\mathbb{R}^{3\times 3}$ using all valid pixels in the
window.

####  III-D2 Pose Optimization

Once the correspondences are established, as shown in Fig. 4, we proceed to
the registration step. We employ the GICP [15] cost function in the
optimization procedure,

$\displaystyle
r^{2}=\|\mu^{\mathcal{P}}-T^{\mathcal{P}}_{\mathcal{S}}\mu^{\mathcal{S}}\|^{2}_{\Sigma},\quad\Sigma=\Sigma^{\mathcal{P}}+R^{\mathcal{P}}_{\mathcal{S}}\Sigma^{\mathcal{S}}R^{\mathcal{S}}_{\mathcal{P}}$
(4)

where $\|\cdot\|_{\Sigma}$ is the Mahalanobis distance and
$\Sigma^{\mathcal{P}}$ and $\Sigma^{\mathcal{S}}$ are the covariance matrices
in the sweep and pano frame respectively. When IMU data is available, we add
an IMU preintegration factor [27] to constrain the start and the end state of
the local trajectory. The error state parameters we wish to optimize represent
a small correction $\Delta T^{\mathcal{P^{\prime}}}_{\mathcal{P}}$ in the pano
frame such that the new state $X^{\mathcal{P^{\prime}}}_{\mathcal{S}}$ can be
updated via

$\displaystyle T^{\mathcal{P^{\prime}}}_{\mathcal{S}}=\Delta
T^{\mathcal{P^{\prime}}}_{\mathcal{P}}\cdot
T^{\mathcal{P}}_{\mathcal{S}},\quad
v^{\mathcal{P^{\prime}}}_{\mathcal{S}}=v^{\mathcal{P}}_{\mathcal{S}}+\Delta
p^{\mathcal{P^{\prime}}}_{\mathcal{P}}/\Delta t$ (5)

where $\Delta p^{\mathcal{P^{\prime}}}_{\mathcal{P}}$ is the translation part
of $\Delta T^{\mathcal{P^{\prime}}}_{\mathcal{P}}$ and $\Delta t=t1-t0$. After
a few iterations of the Levenberg-Marquardt procedure, the resulting update is
applied to the starting state of the sweep trajectory. We then repropagate the
entire trajectory and update correspondences for the next round of ICP,
similar to VICP [28].

###  III-E Map Update

####  III-E1 Depth Fusion

Once the ICP procedure has converged, we can undistort the current sweep using
the optimized local trajectory and update the pano. This is achieved by
projecting every pixel from the sweep to the pano using Eq. 3. Although in
principle one could maintain a probabilistic depth filter for each pixel and
perform recursive Bayesian updates [29] given the new depth, we found that in
practice, the similarity filter used in [22] worked exceptionally well with
very low computation and storage overhead. We keep a counter $k$ alongside the
depth value for each pixel in the pano. If the new depth $d^{\prime}$ is close
enough to the stored depth $d$, we do a weighted average update and increment
the counter by 1, saturating at $k_{\mathrm{max}}$.

$\displaystyle d\leftarrow\frac{k\cdot d+d^{\prime}}{k+1},\quad
k\leftarrow\min(k_{\mathrm{max}},k+1)$ (6)

We otherwise treat the new depth as an outlier and simply decrements the
counter by 1 without modifying the current value. If $k$ reaches 0, we replace
$d$ with $d^{\prime}$. This simple depth fusion strategy produces smooth depth
surfaces as is evident in Fig. 5. It also makes our map resilient to ephemeral
objects, while still being able to adapt to long-term environment changes. For
example, a car driving across the lidar will typically leave a trail of ghost
points in a point cloud map, but it will not be registered into the pano
unless it has stopped moving for a long enough time.

Figure 5: An indoor space (top) mapped with our lidar odometry. The points are
generated by back-projecting a depth panorama. Details of the scene (bottom)
are faithfully preserved by the depth panorama.

####  III-E2 Rendering New Panoramas

As the sensor moves away from the pano, the number of matches between the two
will gradually decrease. This is due to the projective nature of range images,
where data association is affected by occlusions and non-overlapping FOVs.
Therefore, we need to move the pano to a new position in order to retain a
healthy number of matches for ICP.

The decision on whether to relocate the pano is made by looking at the
matching quality $Q$, which is defined as the ratio between the number of
final matches and matching candidates. This number essentially reflects the
amount of non-occluding overlap between the sweep and the pano. Once $Q$ drops
below a certain threshold (0.9 by default), we will render a new pano at the
current location. If IMU measurements are available and motion in roll and
pitch is limited, the new pano will be placed in a gravity aligned frame.
Since we do not have dense depth information in the new location, we do a
forward warping from the old pano to the new one. This will naturally result
in many holes in the new pano, but this can be alleviated by immediately
incorporating the current sweep. The rendering process is executed at a much
lower rate, since it is only triggered when the two viewpoints start to
deviate from each other.

##  IV Implementation

Figure 6: Illustration of how partial sweeps are handled. The incoming partial
sweep is pushed into a sweep buffer. The part being ejected from the buffer is
added to the pano.

###  IV-A Processing Partial Sweeps

To handle partial sweeps, we maintain a circular buffer for a full sweep, and
push incoming partial sweeps into the buffer as they arrive. The part of the
buffer being ejected is immediately fused with the pano using the current
state estimate. This is shown in Fig. 6, where the ejected part is outlined in
red. However, performing ICP with a partial sweep is challenging because the
smaller the size, the more likely it is to encounter degeneracy. Therefore, we
use the entire sweep in the buffer for this task, which has 360° FOV. This
adaptation is key to minimizing latency. For example, if we process half of a
sweep at a time, this reduces the sensor integration time to
$50\text{\,}\mathrm{m}\mathrm{s}$ (assuming a
$10\text{\,}\mathrm{H}\mathrm{z}$ spinning rate). It also halves the feature
detection and pano update time, as they only need to process half of the data.
The ICP time remains unchanged, since it always tries to align the full sweep
in the buffer to the pano. This essentially increases the odometry frequency
to $20\text{\,}\mathrm{H}\mathrm{z}$. One could further divide the sweep into
smaller chunks, but there exists a lower bound on how narrow a partial sweep
can be handled without dropping data, which is subject to the ICP processing
time. In practice, our system can handle 1/8th of a sweep which produces
odometry at $80\text{\,}\mathrm{H}\mathrm{z}$.

###  IV-B Memory Layout and Parallelization

Our implementation follows a data-oriented design [30] and is optimized for
cache-locality and parallelism. We only use array-based data structures and
pre-allocate all storage ahead of time. There is no dynamic allocation
occurring at runtime within our system, which makes it favorable for embedded
applications. Specifically, we only need to allocate enough storage for a full
sweep, a feature grid, two panos and various arrays for the solver and the
local trajectory. Together they require no more than 5-10MB of storage.

The projective data association is also significantly faster compared to
systems that leverage even a carefully designed KD-tree [31]. A tree data
structure requires frequent node allocation and pointer chasing during search
and insertion. Without advanced memory management or allocator support, almost
every node access is a cache miss, which severely hinders performance on
modern processors. On the other hand, the projection operation is constant
time, which is conceptually similar to a hash-table lookup.

Finally, almost every step of our system can be parallelized, since we mostly
operate on images. We choose a single row in both the sweep and the pano as
our atomic work unit, and use a task-based threading library [32] for
scheduling. A finer grain size of a single pixel is possible, but it requires
additional synchronization overhead to avoid false-sharing. Rendering a new
pano is carried out in the background via double buffering
($<$$10\text{\,}\mathrm{m}\mathrm{s}$ in a single thread on Intel CPU), which
does not interfere with the odometry.

##  V Results

We evaluate both the accuracy and runtime performance of our proposed system.
For comparison we choose Fast-LIO2 [6], which is a state-of-the-art lidar-
inertial odometry with an open-source implementation111https://github.com/hku-
mars/FAST_LIO. We use their default configuration for Ouster lidar
222https://ouster.com/products/os1-lidar-sensor with slightly tuned imu noise
parameters. For our system, unless otherwise specified, we use a pano of size
$256\times 1024$, and a cell size of $2\times 16$.

###  V-A Odometry Accuracy

(a) Fast-LIO2

(b) LLOL

Figure 7: Trajectory error plots of Fast-LIO2 and LLOL on Newer College dataset sequence 1. TABLE I: Translational APE and RPE on NC01. Newer College Sequence 1 | Fast-LIO2 | LLOL
---|---|---
Segment | $t_{\mathrm{start}}$ | $t_{\mathrm{end}}$ | APE | RPE | APE | RPE
1 | 25 | 340 | 0.363 | 0.214 | 0.396 | 0.283
2 | 340 | 705 | 0.165 | 0.112 | 0.499 | 0.214
3 | 705 | 1115 | 0.277 | 0.171 | 0.706 | 0.250
4 | 1115 | 1500 | 0.311 | 0.121 | 0.912 | 0.238
TABLE II: Translational APE and RPE on NC05. Newer College Sequence 5 | Fast-LIO2 | LLOL
---|---|---
Segment | $t_{\mathrm{start}}$ | $t_{\mathrm{end}}$ | APE | RPE | APE | RPE
1 | 25 | 100 | 0.104 | 0.172 | 0.130 | 0.328
2 | 100 | 200 | 0.136 | 0.206 | 0.183 | 0.319
3 | 200 | 300 | 0.102 | 0.122 | 0.178 | 0.258
4 | 300 | 400 | 0.162 | 0.119 | 0.262 | 0.288

####  V-A1 Newer College Dataset

As our system consumes lidar packets, this limits our choice of publicly
available datasets to evaluate on. The Newer College dataset [33] provides raw
lidar packets along with ground truth trajectories. We evaluate on two
sequences: 01-short-experiment (NC01) and 05-quad-with-dynamics (NC05), both
collected with a handheld sensor platform at normal walking pace. Note that
NC05 features 4 loops with increasingly aggressive motion. Within each
sequence, the sensor is stopped at various locations for a short while. This
creates well-defined segments which allows for sub-trajectory evaluation. We
report results on these segments in Table I and II, and visualize the full
trajectory of NC01 in Fig. 7.

We report the translational part of Absolute Pose Error (APE) and Relative
Pose Error (RPE) metrics, both generated using the evo [34] package. APE
measures the difference in pose between the reference and the estimated
trajectory. We use Umeyama alignment [35] as a pre-processing step to align
the trajectories for all methods. RPE compares pairwise relative pose within a
small window ($1\text{\,}\mathrm{m}$), which can be used to quantify odometry
drift.

Fast-LIO2 outperforms LLOL in both sequences. This is because part of the
Newer College dataset was collected in a park with many trees and few man-made
structures. This results in reduced planar features for our system. In
addition, the sensor went under tree canopies several times, the heavy
occlusion triggered more frequent rendering of the pano, thus causing extra
drift.

####  V-A2 Our Own Dataset

(a) Seq 1: small loop

(b) Seq 2: large loop

(c) Seq 3: open space

(d) Seq 4: aggressive motion

Figure 8: Qualitative results of Fast-LIO2 (top) and LLOL (bottom) on our own datasets. TABLE III: Trajectory drift of Fast-LIO2 and LLOL on our own dataset. | Fast-LIO2 | LLOL
---|---|---
Seq. | Distance | Error | Drift | Distance | Error | Drift
1 | 69.53 | 0.91 | 1.31% | 64.88 | 0.38 | 0.58%
2 | 257.52 | 6.53 | 2.54% | 244.00 | 0.25 | 0.10%
3 | 208.20 | 0.21 | 0.10% | 201.00 | 0.52 | 0.26%
4 | 869.87 | 358.1 | 41.2% | 63.07 | 0.43 | 0.68%

In addition to the Newer College dataset, we also collect our own for
evaluation. We acquire the dataset with a handheld Ouster OS1 Gen2 sensor
running at $10\text{\,}\mathrm{H}\mathrm{z}$ with 1024 firings per sweep.
Following [36], we come back to the starting point after a loop, and the
difference between the first and the last pose can be used to evaluate
odometry drift without ground truth. We test in four distinct scenarios: a
small loop, a large loop, an open space trajectory and an aggressive motion
trajectory. Quantitative results are shown in Table III. We report the total
distance estimated by each system and their final positional error. The drift
is calculated as the quotient of the two. Both systems work well on the first
3 sequences with normal walking motion, but Fast-LIO2 failed to converge on
the last one with aggressive rotation. LLOL is able to achieve $<$1% drift on
all sequences. Qualitative results are shown in Fig. 8.

###  V-B Runtime and Latency

We benchmark both LLOL and Fast-LIO2 on a desktop CPU (Intel<EMAIL_ADDRESS>as well as an embedded one<EMAIL_ADDRESS>NVIDIA Xavier) under single and
multi-threaded settings. In addition, we conducted experiments on both indoor
and outdoor datasets. The outdoor environment is roughly twice the size of the
indoor one. This is to demonstrate the effect of workspace sizes on each
system. For Fast-LIO2, we report the total time in the same way as the
original paper. For LLOL, we record time from the (partial) sweep arrival to
the end of pano update. The pano rendering only runs infrequently in the
background and is not included in the timing results.

Runtime on our own datasets are shown in Table IV. We see that LLOL is in
general much faster than Fast-LIO2, with sub-millisecond processing time in
some cases. Fast-LIO2 exhibits large variance in runtime between indoor and
outdoor environments, while LLOL runs in near constant time regardless of
workspace size. We also visualize the per-frame runtime on NC01 in Fig. 9.
While LLOL maintains a fairly stable performance over the entire sequence,
Fast-LIO2 displays inconsistent results, with intermittent bursts up to 100ms.
The latency of each method can be computed by adding the sensor integration
time with the runtime. For example, running LLOL with a quarter sweep will
result in an average latency of $28.37\text{\,}\mathrm{m}\mathrm{s}$ on Xavier
with a single thread, as opposed to $189.81\text{\,}\mathrm{m}\mathrm{s}$ for
Fast-LIO2.

Finally, ours system running on embedded CPUs compares favourably to UPSLAM
[22] ($5\text{\,}\mathrm{m}\mathrm{s}$ on Xavier) which requires a GPU.
LITAMIN2 [16] ($5\text{\,}\mathrm{m}\mathrm{s}$ on Intel CPU) can also achieve
similar speed as ours, but the lack of open source implementations precludes a
fair comparison.

TABLE IV: Runtime comparison of Fast-LIO2 vs LLOL on different types of processors under single and multi-threaded settings for both indoor (top) and outdoor (bottom) environments. Results are averaged within the entire sequence. Runtime [ms] | Xavier (ARMv8.2) | Intel (i7-6700K)
---|---|---
method | sweep (ms) | single | multi (8) | single | multi (8)
Fast-LIO2 | full (100) | 44.59 | 27.13 | 22.22 | 21.09
LLOL | full (100) | 12.56 | 4.46 | 5.24 | 1.93
LLOL |  $1/2$ (50) | 6.87 | 3.02 | 3.02 | 1.29
LLOL |  $1/4$ (25) | 4.32 | 2.30 | 1.87 | 0.94
Fast-LIO2 | full (100) | 89.81 | 34.40 | 42.67 | 28.63
LLOL | full (100) | 10.15 | 3.77 | 3.89 | 1.56
LLOL |  $1/2$ (50) | 5.68 | 2.55 | 2.33 | 1.07
LLOL |  $1/4$ (25) | 3.37 | 1.89 | 1.46 | 0.77
Figure 9: Per-frame runtime of Fast-LIO2 and LLOL on NC01, measured on Intel
CPU with 8 threads. LLOL has an average runtime of
1.70$\pm$$0.55\text{\,}\mathrm{m}\mathrm{s}$. Fast-LIO2 has an average runtime
of 21.73$\pm$$16.50\text{\,}\mathrm{m}\mathrm{s}$.

##  VI Conclusions

In this paper, we present a low-latency odometry system for spinning lidars.
It has comparable accuracy and drift to a state-of-the-art lidar odometry
system [6], but runs an order of magnitude faster. The decision to handle
partial sweeps significantly boosts the odometry update frequency of our
system compared to approaches operating at the nominal rotational rate of the
lidar sensor. A key enabler of our fast implementation is the use of depth
panoramas as our local map representation, this affords us constant time
search and update complexity across different environment scales. Our odometry
system offers a high data throughput combined with a low memory footprint,
these features make it suitable for deployment on resource constrained
platforms. We open source our implementation as a contribution to the
community.

## References

  * [1] J. Zhang and S. Singh, “Visual-lidar odometry and mapping: low-drift, robust, and fast,” _2015 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 2174–2181, 2015.
  * [2] T. Shan, B. Englot, C. Ratti, and D. Rus, “Lvi-sam: Tightly-coupled lidar-visual-inertial odometry via smoothing and mapping,” _ArXiv_ , vol. abs/2104.10831, 2021.
  * [3] S. Zhao, H. Zhang, P. Wang, L. Nogueira, and S. Scherer, “Super odometry: Imu-centric lidar-visual-inertial estimator for challenging environments,” _ArXiv_ , vol. abs/2104.14938, 2021.
  * [4] D. Frossard, S. Suo, S. Casas, J. Tu, R. Hu, and R. Urtasun, “Strobe: Streaming object detection from lidar packets,” _ArXiv_ , vol. abs/2011.06425, 2020.
  * [5] W. Han, Z. Zhang, B. Caine, B. Yang, C. Sprunk, O. Alsharif, J. Ngiam, V. Vasudevan, J. Shlens, and Z. Chen, “Streaming object detection for 3-d point clouds,” _ArXiv_ , vol. abs/2005.01864, 2020.
  * [6] W. Xu, Y. Cai, D. He, J. Lin, and F. Zhang, “Fast-lio2: Fast direct lidar-inertial odometry,” _ArXiv_ , vol. abs/2107.06829, 2021.
  * [7] J. Zhang and S. Singh, “Loam: Lidar odometry and mapping in real-time,” in _Robotics: Science and Systems_ , 2014.
  * [8] T. Shan and B. Englot, “Lego-loam: Lightweight and ground-optimized lidar odometry and mapping on variable terrain,” _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 4758–4765, 2018.
  * [9] Z. Liu and F. Zhang, “Balm: Bundle adjustment for lidar mapping,” _IEEE Robotics and Automation Letters_ , vol. 6, pp. 3184–3191, 2021.
  * [10] C. Park, P. Moghadam, S. Kim, A. Elfes, C. Fookes, and S. Sridharan, “Elastic lidar fusion: Dense map-centric continuous-time slam,” _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 1206–1213, 2018\.
  * [11] C. Park, P. Moghadam, J. Williams, S. Kim, S. Sridharan, and C. Fookes, “Elasticity meets continuous-time: Map-centric dense 3d lidar slam,” _ArXiv_ , vol. abs/2008.02274, 2020.
  * [12] J. Behley and C. Stachniss, “Efficient surfel-based slam using 3d laser range data in urban environments,” in _Robotics: Science and Systems_ , 2018.
  * [13] J. Quenzel and S. Behnke, “Real-time multi-adaptive-resolution-surfel 6d lidar odometry using continuous-time trajectory optimization,” _ArXiv_ , vol. abs/2105.02010, 2021.
  * [14] M. Yokozuka, K. Koide, S. Oishi, and A. Banno, “Litamin: Lidar-based tracking and mapping by stabilized icp for geometry approximation with normal distributions,” _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 5143–5150, 2020.
  * [15] A. V. Segal, D. Hähnel, and S. Thrun, “Generalized-icp,” in _Robotics: Science and Systems_ , 2009.
  * [16] M. Yokozuka, K. Koide, S. Oishi, and A. Banno, “Litamin2: Ultra light lidar-based slam using geometric approximation applied with kl-divergence,” _ArXiv_ , vol. abs/2103.00784, 2021.
  * [17] K. Li, M. Li, and U. Hanebeck, “Towards high-performance solid-state-lidar-inertial odometry and mapping,” _IEEE Robotics and Automation Letters_ , vol. 6, pp. 5167–5174, 2021.
  * [18] J. Lin and F. Zhang, “Loam livox: A fast, robust, high-precision lidar odometry and mapping package for lidars of small fov,” _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 3126–3131, 2020\.
  * [19] T. Shan, B. Englot, D. Meyers, W. Wang, C. Ratti, and D. Rus, “Lio-sam: Tightly-coupled lidar inertial odometry via smoothing and mapping,” _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 5135–5142, 2020.
  * [20] H. Ye, Y. Chen, and M. Liu, “Tightly coupled 3d lidar inertial odometry and mapping,” _2019 International Conference on Robotics and Automation (ICRA)_ , pp. 3144–3150, 2019.
  * [21] W. Xu and F. Zhang, “Fast-lio: A fast, robust lidar-inertial odometry package by tightly-coupled iterated kalman filter,” _IEEE Robotics and Automation Letters_ , vol. 6, pp. 3317–3324, 2021.
  * [22] A. Cowley, I. D. Miller, and C. J. Taylor, “Upslam: Union of panoramas slam,” _ArXiv_ , vol. abs/2101.00585, 2021.
  * [23] C. J. Taylor, A. Cowley, R. Kettler, K. Ninomiya, M. Gupta, and B. Niu, “Mapping with depth panoramas,” _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 6265–6272, 2015.
  * [24] G. Kim and A. Kim, “Scan context: Egocentric spatial descriptor for place recognition within 3d point cloud map,” _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 4802–4809, 2018.
  * [25] H. Wang, C. Wang, and L. Xie, “Intensity scan context: Coding intensity and geometry relations for loop closure detection,” _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 2095–2101, 2020\.
  * [26] J. Serafin and G. Grisetti, “Nicp: Dense normal based point cloud registration,” _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 742–749, 2015.
  * [27] C. Forster, L. Carlone, F. Dellaert, and D. Scaramuzza, “On-manifold preintegration for real-time visual–inertial odometry,” _IEEE Transactions on Robotics_ , vol. 33, pp. 1–21, 2017.
  * [28] S. Hong, H. Ko, and J. Kim, “Vicp: Velocity updating iterative closest point algorithm,” _2010 IEEE International Conference on Robotics and Automation_ , pp. 1893–1898, 2010.
  * [29] C. Forster, M. Pizzoli, and D. Scaramuzza, “Svo: Fast semi-direct monocular visual odometry,” _2014 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 15–22, 2014.
  * [30] R. Fabian, _Data-Oriented Design: Software Engineering for Limited Resources and Short Schedules_ , 2018.
  * [31] Y. Cai, W. Xu, and F. Zhang, “ikd-tree: An incremental k-d tree for robotic applications,” _ArXiv_ , vol. abs/2102.10808, 2021.
  * [32] M. Voss, R. Asenjo, and J. Reinders, _Pro TBB: C++ Parallel Programming with Threading Building Blocks_ , 1st ed. USA: Apress, 2019.
  * [33] M. Ramezani, Y. Wang, M. Camurri, D. Wisth, M. Mattamala, and M. Fallon, “The newer college dataset: Handheld lidar, inertial and vision with ground truth,” _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 4353–4360, 2020.
  * [34] M. Grupp, “evo: Python package for the evaluation of odometry and slam.” https://github.com/MichaelGrupp/evo, 2017.
  * [35] S. Umeyama, “Least-squares estimation of transformation parameters between two point patterns,” _IEEE Transactions on Pattern Analysis & Machine Intelligence_, vol. 13, no. 04, pp. 376–380, 1991.
  * [36] J. Engel, V. Usenko, and D. Cremers, “A photometrically calibrated benchmark for monocular visual odometry,” _ArXiv_ , vol. abs/1607.02555, 2016.