# Cold atoms in low dimensions – a laboratory for quantum dynamics

S. I. Mistakidis  A. G. Volosniev  R. E. Barfknecht  T. Fogarty  Th. Busch  A.
Foerster  P. Schmelcher  N. T. Zinner  ITAMP, Center for Astrophysics $|$
Harvard $\&$ Smithsonian, Cambridge, MA 02138 US  IST Austria (Institute of
Science and Technology Austria), Am Campus 1, 3400 Klosterneuburg, Austria
Instituto de Física da UFRGS, Av. Bento Gonçalves 9500, Porto Alegre, RS,
Brazil  Quantum Systems Unit, Okinawa Institute of Science and Technology
Graduate University, Okinawa, Japan 904-0495  Center for Optical Quantum
Technologies, University of Hamburg, Department of Physics, Luruper Chaussee
149, D-22761, Hamburg, Germany  The Hamburg Centre for Ultrafast Imaging,
University of Hamburg, Luruper Chaussee 149, D-22761, Hamburg, Germany
Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120,
DK-8000 Aarhus C, Denmark

###### Abstract

Cold atomic gases have become a paradigmatic system for exploring fundamental
physics, while at the same time allowing to explore applications in quantum
technologies. The accelerating developments in the field have led to a highly
advanced set of engineering techniques that, among others, can tune
interactions, shape the external geometry, select among a large set of atomic
species with different properties, as well as control the number of atoms. In
particular, it is possible to operate in lower dimensions and drive atomic
systems in the strongly correlated regime. In this review, we discuss recent
advances in few-body cold atom systems confined in low dimensions from a
theoretical viewpoint. We mainly focus on bosonic systems in one dimension and
provide an introduction to the well-known static properties before we discuss
the state-of-the-art quantum dynamical processes stimulated by the presence of
correlations. Besides describing the fundamental physical phenomena arising in
these systems, we also provide an overview of the tools and methods that are
commonly used, thus delivering a balanced and comprehensive overview of the
field. We conclude by giving an outlook on possible future directions that are
interesting to explore in these correlated systems.

###### keywords:

Cold atoms, one-dimensional systems, few-body physics, multicomponent bosonic
settings, quantum dynamics, droplets, spinor gases, systems with impurities,
quantum control, long-range physics, beyond mean-field approaches

††journal: Physics Reports

###### Contents

  1. 1 Introduction
    1. 1.1 Cold atom platforms
    2. 1.2 Reducing the dimension
    3. 1.3 Towards the few-body limit
    4. 1.4 Tuning the interaction
    5. 1.5 Review outline
  2. 2 Fundamental models and exact solutions
    1. 2.1 Bose gas with contact interactions
    2. 2.2 Solutions in the few-body limit
    3. 2.3 The Tonks-Girardeau model
    4. 2.4 Fermionic systems
  3. 3 Single component static and dynamics
    1. 3.1 Two-particles in a harmonic trap
    2. 3.2 Three particles
    3. 3.3 Approximate approaches
    4. 3.4 Few-particle properties
    5. 3.5 Dynamical properties and quenches; time-dependent modulations
  4. 4 Two-component 1D gases
    1. 4.1 Ground state properties of bosonic mixtures
    2. 4.2 Energy transfer and breathing modes
    3. 4.3 Quench dynamics and pattern formation
    4. 4.4 Tunneling dynamics
    5. 4.5 Tunneling with mass imbalance
    6. 4.6 Bose-Fermi mixtures
    7. 4.7 Systems close to fermionization
  5. 5 Impurities and polarons
    1. 5.1 A single impurity
    2. 5.2 Two impurities
    3. 5.3 Further directions
  6. 6 Spin-1 bosonic gases
    1. 6.1 Overview and state-of-the-art
    2. 6.2 One-dimensional setting and basic definitions
    3. 6.3 Mean-field picture
    4. 6.4 Phase diagram
    5. 6.5 Spin-mixing dynamics and the single mode approximation
  7. 7  Control: driven dynamics and transport
    1. 7.1 Driven dynamics
    2. 7.2 Transport
    3. 7.3 Quantum control
  8. 8 Quantum droplets
    1. 8.1 The generalized Gross-Pitaevskii framework
    2. 8.2 Stationary properties of one-dimensional droplets
    3. 8.3 Dynamics of one-dimensional droplets
  9. 9 Methods: review
    1. 9.1 Exact Diagonalization
    2. 9.2 Bethe Ansatz-Based Methods
    3. 9.3 Density Matrix Renormalization Group
    4. 9.4 The Multi-Layer Multi-Configuration Time-Dependent Hartree Method for Atomic Mixtures
    5. 9.5 Spin-Chain methods
  10. 10 Beyond contact interactions
    1. 10.1 Finite-range corrections in atom-atom interactions
    2. 10.2 Odd-wave interactions
    3. 10.3 Dipole-dipole interactions
    4. 10.4 Dipolar droplets
    5. 10.5 Atom-ion interactions
  11. 11 Summary and future prospectives

##  1 Introduction

Cold atom experiments have over the last two decades emerged as a very strong
contender for the study of fundamental physics and for the simulation of
quantum behavior relevant for other branches of physics (Bloch et al., 2008;
Zohar et al., 2015; Saffman, 2016; Gross and Bloch, 2017; Szigeti et al.,
2021). One particular intriguing problem that has been addressed is the
quantum many-body problem, which has excited researchers and students over
decades. It is, however, a distinctly difficult problem, as the complex
dynamics such systems show is governed by an exponentially large number of
degrees of freedom that can rarely be fully accounted for. Nevertheless over
the past century many successful approaches were developed, many of which have
become textbook material, and intriguing physics has been discovered.

While traditionally solid state systems have been the driver of quantum many-
body physics, trapped ultracold atoms have recently emerged as a paradigmatic
system to study fundamental questions of this area. This is due to the large
amount of experimental control one can have over these systems as well as
their purity, which often makes them amenable to accurate theoretical
calculation. In particular, the experimental realization of small systems in
lower dimensions, that allow for strong correlations and an extensive control
over almost all degrees of freedom, boosted this field. This progress has made
it feasible to study the crossover from few- to many-body physics, providing a
route to study quantum-many body systems starting from few-body correlations
in a bottom-up approach.

In addition to their value for contributing to the understanding of the
quantum many-body problem, small samples have been recognised as promising
systems by themselves, in which, for example, one can quantum engineer
specific non-classical physics. The high level of control available in cold
atom experiments can be used to explore a diverse range of few-body quantum
effects, including the generation and control of entanglement, spin-dependent
transport, the dynamical orthogonality catastrophe, few-body quantum ‘chaos’,
filamentation processes, dynamical phase separation, interaction-induced
tunneling and Josephson-type dynamics, spontaneous droplet formation and many
more.

In this review we summarise the recent body of work on the theory of few-body
dynamics and the connections to many-body behavior in cold gases. Our emphasis
is on setups and situations which are addressable in present and near-future
experiments, thereby hoping that our discussion may stimulate further
experimental and theoretical efforts. It is worth noting that other useful
reviews exist on several topics related to the use of cold atoms as quantum
simulators for many-body condensed-matter systems (Cazalilla et al., 2011;
Guan et al., 2013). The reviews by Blume (2012) and by Sowiński and García-
March (2019) provide good introductions to the physics of few-body systems in
low-dimensional cold-atom setups, but have limited overlap with what will be
covered in the present review, as they mainly focus on static phenomena.

While we focus on theoretical work and approaches, we will in the following
also briefly account for selected aspects of the experimental progress made in
cold atom systems in recent years, as this is closely related and often the
motivation for many of the theoretical developments that we will review in the
following sections.

###  1.1 Cold atom platforms

One advantage of having tunable trapping potentials is the possiblity to
freeze out dynamics in certain spatial regions through strong confinement.
While the area of cold gases was initially limited to weakly interacting
systems, a main driver of low-dimensional setups has been the possibility to
go to different correlation regimes. This brought the field into a position
where some famous model predictions could be experimentally realized and
tested.

Arguably, one of the main drivers behind the early pursuits of low-dimensional
setups was the seminal results of Girardeau (Girardeau, 1960), who pointed out
that strongly interacting Bose systems in one dimension display fermionic
aspects in their behavior. In the early 2000s, this so-called Tonks-Girardeau
regime became a topic of great interest due to the ability to produce
controllable samples of cold atoms in lower-dimensional geometries (Bloch et
al., 2008). For the theoretical treatment of these systems, a key step was the
work by Olshanii (Olshanii, 1998) that provided a comprehensive theoretical
model of how to simulate 1D physics, including very strong interactions, by
tuning the confining geometry of the system. At the same time the first
experiments capable of tuning the atom-atom interactions in cold gas samples
were reported (Inouye et al., 1998), and tunability using the so-called
Feshbach resonance technique (Chin et al., 2010) is now part of the standard
toolbox of cold atom experiments. Combining the tunability of the geometry and
the interaction strength results in so-called confinement-induced resonances
111The physical origin of the CIRs is that they fulfill a Fano-Feshbach
mechanism where all the energetically closed transversal modes of the trapping
potential collectively support one bound state that lies in the energy
interval of the open channel (Friedrich and Friedrich, 2006; Greene et al.,
2017), i.e. the ground state transverse mode.  (CIR) (Bergeman et al., 2003).
These were theoretically predicted to provide an essentially zero-range
interaction with a coupling strength that can be tuned to large positive or
negative values through the adjustment of the magnetic field used in Feshbach
resonances. This was a key feature of low-dimensional geometries and implied
that cold atomic gases provide a universal platform for studying some of the
many one-dimensional models that have been proposed in other branches of
physics (Giamarchi, 2003; Sutherland, 2004) in a very direct and tunable
manner. Hence, it was no surprise that several experimental groups set out to
test these predictions very soon after.

###  1.2 Reducing the dimension

The physics that can be observed in reduced dimensions has a number of aspects
that depend on whether one looks at few- or many-body systems. The key to
understanding these different behaviours is a comparison of the kinetic and
interaction energies as function of spatial dimension. A simple picture to use
here is that of a particle in a box or in a harmonic trap. As known from
elementary quantum mechanics, the zero-point energy increases as the number of
spatial dimensions increases. On the one hand, this implies that to correlate
few-body systems one must use stronger interactions in 3D as compared to 2D or
1D, and indeed this can be shown by elementary methods to be the case (Landau
and Lifshitz, 1977). On the other hand, for larger systems, i.e. going towards
the many-body regime, different dynamics appear, and certain collective
behavior in large systems is actually harder to achieve in 1D and 2D as
compared to 3D. An example of this is Bose-Einstein condensation, which is
considerably harder to achieve (and even define) in 1D and 2D (Pethick and
Smith, 2002). The foundation of this result is the famous Mermin-Wagner
theorem that excludes long-range order in lower dimensions (Mermin and Wagner,
1966), and hence prohibits many forms of collective behavior. Note that this
is strictly the case in homogeneous systems only, and the presence of traps
can lead to, for instance, quasi-condensates even in 2D and 1D (Pethick and
Smith, 2002). These theoretical arguments imply that changes in dimensionality
can lead to vastly different physical behavior both in few- and many-body
physics, and are a strong motivation for experimental work in cold atoms with
low-dimensional geometries.

Figure 1:  Measurement of a 3D to 1D crossover through the excitation spectrum
of an array of 1D Bose gases as function of the depth of the trapping
potential quantified through the ratio of interaction to hopping, $U/J$. The
appearance of a discrete structure due to a vanishing superfluid component is
clearly visible. The 3D case shows the superfluid to Mott insulator
transition. Image from (Stöferle et al., 2004)

Some of the first experiments to approach a one-dimensional geometry were
using a combination of optical and magnetic traps, such that the motion of
atoms along the transverse directions in a weakly-interacting Bose-Einstein
condensate were restricted and effectively frozen out (Görlitz et al., 2001;
Dettmer et al., 2001). Alternatively, it was also shown how a combined system
of fermions and bosons can arrange itself in a magnetic trap so as to achieve
an effective one-dimensional confinement of the Bose gas component (Schreck et
al., 2001). Other ways to achieve geometric confinement used optical lattices
(Greiner et al., 2001; Moritz et al., 2003; Tolra et al., 2004) but only
intermediate confinement strengths could be achieved. This changed drastically
in 2004, when two seminal experiments (Paredes et al., 2004; Kinoshita et al.,
2004) managed to reach the limit of very strong interactions, the so-called
Tonks-Girardeau regime, using a novel approach in which a (shallow) optical
lattice was applied along the 1D spatial direction. The shallow longitudinal
lattice quenches the kinetic energy of the particles but not enough for them
to become trapped in the minima. Rather, the effect is to increase the
effective mass of the bosons in the experiments, and this means that the
relative contribution of the interactions between the particles can dominate
the kinetic term. This has made it possible to study the transition from 3D
down to 1D physics with strong interactions (Stöferle et al., 2004) (see Fig.
1). These experiments provided the verification of the fact that one-
dimensional Bose gases do indeed adopt a fermionic nature in the presence of
strongly (repulsive) short-range interactions as predicted by (Girardeau,
1960; Lieb and Liniger, 1963). They also provided a stepping stone to cold-
atom quantum simulators of low-dimensional many-body models. Moreover, it was
possible to probe different features including the wave function overlap in
the system (Kinoshita et al., 2005), and also the non-equilibrium dynamics in
the Tonks-Girardeau regime (Kinoshita et al., 2006).

The 1D Bose gas is a topic of continued pursuit and in recent years the so-
called Tomonaga-Luttinger-Liquid behavior in quantum critical systems has been
observed in cold atoms (Yang et al., 2017) (similar physics was also observed
in a superconducting Josephson array (Cedergren et al., 2017)), which is a
prediction stemming from condensed-matter theory. Furthermore, a recent study
by (Wilson et al., 2020) showed how the 1D Bose gas in the strongly
interacting regime will show dynamical fermionization in its momentum
distribution as it is released from the trap (Rigol and Muramatsu, 2005;
Minguzzi and Gangardt, 2005). The observation of transport properties in 1D
strongly interacting Bose gases is also a topic of great importance, and it
can be done by injecting an impurity in the gas and studying its motion
(Palzer et al., 2009; Rutherford et al., 2011). A recent experiment along this
line was reported by (Meinert et al., 2017), who showed that with strong
interactions in a 1D Bose gas, a single impurity will display Bloch
oscillations. This effect arises from the background Tonks-Girardeau gas that
acts as an effective lattice for the impurity. Another approach to studying
impurities in the presence of a cold gas of background atoms has been
pioneered by (Spethmann et al., 2012) with experiments that can control
individual impurities within the gas. All of these demonstrate the versatility
of cold bosons as a platform to perform quantum simulation in low-dimensional
geometries, including the successful realization of topological phases in 1D
with Bose gases (Atala et al., 2013; Lohse et al., 2016), and in ladder setups
(Atala et al., 2014; Mancini et al., 2015).

It is worth noting that another direction to achieve low-dimensional
confinement is to use so-called atom chips (Fortagh et al., 1998; Müller et
al., 1999; Reichel et al., 1999; Dekker et al., 2000; Folman et al., 2000).
These are surfaces with carefully crafted nano-wires that can produce a
trapping potential close to the surface for an ensemble of atoms. This has
allowed for a number of very impressive experiments in particular with bosonic
samples and the study of matter-wave interferometry (Schumm et al., 2005; Jo
et al., 2007), 1D models of thermodynamics (van Amerongen et al., 2008), and
quantum and thermal noise (Hofferberth et al., 2007, 2008). More information
about the atom chip developments can be found in dedicated reviews (Fortágh
and Zimmermann, 2007; Folman et al., 2002).

###  1.3 Towards the few-body limit

A seminal experimental achievement was reached around a decade ago when
samples with less than ten particles near zero temperature were created in a
one-dimensional trapping potential with fermionic 6Li atoms (Serwane et al.,
2011). This was done by loading a large number of atoms into a trap with a
small ‘dimple’ at the bottom and then ‘cutting’ the trap to remove most of the
atoms, see Fig. 2. This experimental procedure paved the way for several
studies of the basic properties of quantum systems starting from the few-body
limit and building towards the collective states seen with many particles.
This setup was used to create a simple two-body system in which the
interaction strength could be varied from weak to strong (Zürn et al., 2012).
It also allowed for the explicit experimental demonstration of the predictions
of confinement-induced resonances for the smallest system possible (Olshanii,
1998). The experimental work confirmed that indeed local properties of two
particles with opposite internal (spin) states are identical to those of
spinless fermions in the limit where the effective 1D interaction strength
diverges (Zürn et al., 2012). This was shown by using a novel tunneling
technique in which one particle could ‘spill’ out of the trap with a tunneling
rate that is dependent on the interaction strength of the pair. As a check on
the setup, the experiment has been performed using two fermions in identical
internal states, for which the tunneling rate is independent of the
interaction strength as identical fermions cannot interact via a short-range
$s$-wave interaction (Zürn et al., 2012).

These experimental techniques were later used to study the case of attractive
interactions with more than two particles. It was shown that pairing occurs in
system sizes down to six particles (Zürn et al., 2013) with a clear signature
in terms of an odd-even effect known from pairing studies in nuclei (Ring and
Schuck, 1980; Jensen et al., 2004; Zinner and Jensen, 2013) or quantum dots
(Kouwenhoven et al., 2001; Reimann and Manninen, 2002; Hanson et al., 2007).

Figure 2:  (A) Starting from a degenerate two-component Fermi gas of about 600
atoms in a microtrap, a few-particle system is created by adiabatically
deforming the potential to allow atoms in higher levels to tunnel out. (B) The
mean atom number decreases when the depth of the trap is reduced. This happens
in step of two because each energy level in the trap is occupied with one atom
per spin state. The variance is shown on the right and one can see that for
even atom numbers, the number fluctuations are strongly suppressed. Image from
(Serwane et al., 2011)

The availability of tunable interactions and tunable particle numbers,
prompted an immediate experimental investigation into the emergence of many-
body properties related to the impurity problem. The work of (Wenz et al.,
2013) showed how to study the impurity problem starting from two particles and
building it bottom-up. This effectively corresponded to constructing a Fermi
sea atom by atom, and the energetics could be studied and compared to
theoretical predictions, as will be discussed in more detail later on.
Remarkably this experiment also found that already for the case of five or six
identical fermions interacting with a single particle of opposite spin the
system has an energy in agreement with many-body models for (repulsive)
interaction strengths ranging from weak to strong. More specifically, the 1D
Fermi impurity problem, or polaron problem, has a known solution without a
trapping potential (McGuire, 1965). The experimental results approach the
energy prediction for the polaron in the thermodynamic limit already for only
moderate system sizes, indicating a fast convergence to a full Fermi sea
interaction with an impurity in the 1D system.

Having access to two-component fermions in tailored 1D traps in the strongly
interacting regime enables the study of spin dynamics for few-body systems,
and hence the study of magnetic correlations from the few-body point of view.
In particular, the work of (Murmann et al., 2015a) showed how to use a
tailored trap in 1D to build a Hubbard model starting from just two atoms.
This means that the dynamics of the system as a function of the typical
parameters can be explored. These parameters have typically been taken as
fixed in theoretical Hubbard model studies of the many-body system, and it has
remained uncertain to what extent they represent the true few-body parameters
(Valiente and Petrosyan, 2009; Büchler, 2010; Valiente and Zinner, 2015). By
using a double-well setup with a tunable barrier and a Feshbach resonance to
tune the atom-atom coupling, (Murmann et al., 2015a) could explore the system
going from the few-body equivalent of the Mott insulator to a superfluid
phase, as well as enter the regime with attractive interactions. These kinds
of experiments provide a starting point for building a quantum simulator of
paradigmatic condensed-matter models from the bottom up. Furthermore, in the
case of strong (repulsive) interactions, the Hubbard model becomes effectively
a spin model as the motional degrees of freedom freeze out in 1D. Here the
system should enter an antiferromagnetic regime. The work of (Murmann et al.,
2015b) showed that antiferromagnetic spin correlations could be detected for
strong interactions in systems with three and four particles in a harmonic
trap using the tunneling techniques mentioned above.

Even though our focus here is on small and continuous systems in low
dimensions, it is necessary to briefly mention the Hubbard model and magnetic
correlation studies in optical lattices as they are effectively also probing
few-body correlations. One of the first direct observations of the
antiferromagnetic spin correlations was reported in (Greif et al., 2013) using
an optical lattice setup with fermionic atoms. This was followed by work of
(Boll et al., 2016) on fermionic chains of atoms in optical lattices with
access to both the spin and charge (atomic density) degrees of freedom, and by
work on 2D lattice setups (Mazurenko et al., 2017). 1D Fermi-Hubbard chains of
atoms have since been used to elucidate a lot of interesting physics in
experiments including string ordering away from unit filling (Hilker et al.,
2017), incommensurate spin-density waves (Salomon et al., 2019), signatures of
the dynamical deconfinement of spin and charge (Vijayan et al., 2020) and in
the 2D context the direct observation of magnetic polarons (Koepsell et al.,
2019).

In recent developments, it has become possible to image the low-dimensional
few-body Fermi systems (Bergschneider et al., 2019). This was not only
achieved in the traditional sense of knowing the spatial position of atoms,
but also imaging in momentum space, see Fig. 3. A double-well setup such as
the one in (Murmann et al., 2015a) can then be combined with these new imaging
techniques to study the correlations of as few as two particles and provide
information about the level of entanglement in the system. These imaging
techniques have also enabled the study of controlled few-body collisions in
the system (Guan et al., 2019). This can be used to study details of the
dynamics in collisions and quench dynamics, including the question of how the
system may or may not enter relaxation, opening up a few-body angle on
equilibrium, non-equilibrium and quantum chaotic behavior (Guan et al., 2019).
An extension of the manipulation and imaging techniques to 2D arrays of
fermion pairs was recently reported (Hartke et al., 2022) that provides very
long coherence times for potential use in cold atom-based quantum information
and computing setups.

Figure 3:  (a) Position and (b) momentum correlation functions of a Fermi-
Hubbard dimer trapped in a double well potential for different interaction
strengths, $U/J$, going from the attractive, through non-interacting and onto
the strongly repulsive regime. Image from Bergschneider et al. (2019).

The imaging techniques and control are part of the overall toolbox within the
field of optical tweezers used for cold atoms, and we therefore briefly
discuss a few of the recent results from optical tweezers applied to cold
atoms. The work towards local manipulation of atoms in traps has a long
history (Dumke et al., 2002; Miroshnychenko et al., 2006; Beugnon et al.,
2007; Karski et al., 2009; Grünzweig et al., 2010; Muldoon et al., 2012). In
recent years optical tweezers have seen a tremendous development, in
particular in bosonic samples. The seminal work of (Kaufman et al., 2015)
showed how to use spin-exchange interactions in 87Rb atoms to achieve
position-resolved coherence measurements and entangle two atom systems in a
controlled manner. Furthermore, the same kind of optical tweezer setup can be
used to deterministically create entanglement even in systems where the atoms
are non-interacting (Lester et al., 2018). This was achieved using a beam-
splitter and post-selection technique akin to the Hong-Ou-Mandel effect
(Kaufman et al., 2014; Lopes et al., 2015; Preiss et al., 2015), and can be
seen as a step towards a protocol for measurement-based quantum information
processing using atoms. The few-body results were soon followed by experiments
that showed how large-scale tweezer arrays can be experimentally realized. In
(Endres et al., 2016) a total of 100 optical tweezers was used to obtain 1D
arrays of up to 50 bosonic atoms with no defects in the system in a remarkably
fast protocol, paving the way for structured assembly of 1D atomic systems for
quantum information and quantum simulation purposes. Moreoever, (Endres et
al., 2016) also showed how one may replace lost atoms in the arrays. A similar
feat was accomplished by (Barredo et al., 2016) who used up to 50 microtraps
to assemble an array of bosonic atoms into a 2D lattice with full tunability.
We note that it has even become possible recently to create an array of cold
molecules using an optical tweezer setup (Anderegg et al., 2019) and this may
open up the possibility to have longer-range interactions in such setups.

All of the mentioned tweezer experiments thus far are based on bosonic atoms,
and moving to a fermionic system was a remaining challenge until very
recently. In (Preiss et al., 2019), a system of 1D fermions was captured in a
set of optical tweezers and interference effects with high contrast were
demonstrated. The key measurements involved the momentum distributions that
could be used to reveal two- and three-body correlation functions for small
system sizes (Becher et al., 2020). This allows one to study Hanbury-Brown and
Twiss (anti)-correlations for fermions using cold atomic gases in controlled
setups. Furthermore, in a two-dimensional setting, the ability to probe the
momentum distribution led to the first observation of a so-called ‘Pauli
crystal’ (Holten et al., 2021) . This is an elusive state that features
higher-order density correlations between trapped fermions even without any
interactions. It has also become possible to observe few-body shell effects in
two-dimensional traps with fermions directly (Bayha et al., 2020). For the
two-dimensional setups it has been conjectured that by rotating the trap one
may even access the exotic physics of fractional quantum Hall systems using
either bosonic (Andrade et al., 2021; Palm et al., 2021) or fermionic (Palm et
al., 2020) few-body systems. These new developments show how the imaging and
tweezer techniques can be used also in low-dimensional systems for fermions
and provide very interesting prospects for future work.

###  1.4 Tuning the interaction

Figure 4:  Effective 1D coupling strength $g_{1D}$ (in units of $\hbar/(\mu
a_{\perp})$) with respect to $a_{s}(E)/a_{\perp}$ ($a_{S}(E)$ denotes the
scattering length). The location of the confinement-induced resonance is
indicated by the red star. The solid line depicts the analytical prediction of
Eq. (1), while the black dots refer to numerical calculations of the two boson
problem with van-der-Waals two-body interactions. Image from (Greene et al.,
2017).

One of the key elements behind the success of cold atom experiments is the
possibility to tune interactions using external fields. For quasi-1D systems,
Olshanii (1998) argued that the pairwise interactions can be modeled via a
one-dimensional $\delta$-function interaction whose strength is given by

$g_{1D}=\frac{2\hbar^{2}a_{s}}{\mu a_{\perp}(a_{\perp}-Ca_{s})},$ (1)

where $a_{s}$ is the three-dimensional scattering length, $C$ is a constant,
$a_{\perp}$ denotes the harmonic oscillator length in the direction of the
transverse confinement and $\mu$ is the reduced mass. Equation (1)
demonstrates that by changing $a$, one can adjust both the magnitude and the
sign of the one-dimensional interaction strength. The behavior of the
effective coupling constant $g_{1D}$ as one varies $a_{s}/a_{\perp}$ is shown
in Fig. 4. The analytical result of Eq. (1) (solid orange line) is in
agreement with the numerical solution of a model two-body Hamiltonian with the
potential $V_{vdw}=C_{10}/r^{10}-C_{6}/r^{6}$ (black dots). Notice also that
for $a_{\perp}/a_{s}\to C$ the effective coupling constant diverges,
$g_{1D}\to\infty$. This marks the position of the confinement-induced
resonance (CIR) and is indicated by the red star. A more detailed discussion
on this topic, including also the behavior of the bound states, can be found
in (Greene et al., 2017).

Figure 5:  A: Ratio of the frequency of the lowest axial compression mode over
the lowest dipole mode, R, vs. the interaction strength. The blue squares show
the measurements in the attractive regime and provide evidence for the super-
Tonks-Girardeau gas. The red circles show the transition from the Thomas-Fermi
to the Tonks-Girardeau regime for repulsive interactions. B: The three-
dimensional scattering length $a_{\text{3D}}$ (dashed-dotted), the one-
dimensional scattering length $a_{\text{1D}}$ (solid), and the one-dimensional
interaction parameter $g_{\text{1D}}$ (dashed) in the vicinity of the Feshbach
resonance (FR) at 47.78 G. Image from (Haller et al., 2009)

Tunability of $g_{1D}$ has allowed experimentalists to prepare a one-
dimensional Bose gas in the strongly interacting regime (Haller et al., 2009),
as well as to transfer the gas from the strongly repulsive regime and into a
regime with strong attractive interactions (see Fig. 5). The latter regime has
been denoted the super-Tonks-Girardeau regime (Astrakharchik et al., 2005). It
is important to point out that if one prepares a super-Tonks-Girardeau state
by moving slowly from a regime of strong repulsion into a regime of strong
attraction by crossing the CIR, the resulting state is a highly excited state.
Nevertheless, this state was shown to live long enough for it to be studied
(Haller et al., 2009), and the nature of the CIRs could subsequently be
explored in greater detail (Haller et al., 2010).

Before closing our discussion on interactions, we briefly comment on the
physics associated with fermions atoms in the presence of tunable interactions
in one dimension. Here a clear signature of CIR was already seen in early
experiments on two-fermion molecule formation (Moritz et al., 2005). Probing
the binding energy of these molecules using radio-frequency (RF) spectroscopy
showed that the 1D situation was very different from a 3D one, and that the
overall CIR predictions are applicable to Fermi systems as well. More
recently, the shift in position of the resonances due to the confinement has
even allowed for more accurate determination of the position of the 3D
Feshbach resonance using small samples of fermionic molecules (Zürn et al.,
2013). The corrections due to anharmonicity of the confinement has also been
experimentally addressed using a sample of just two fermions and considering
molecule formation using RF spectroscopy (Sala et al., 2013). The discussion
thus far has focused on $s$-wave Feshbach resonances, but higher partial waves
are also modified due to the confinement (Hess et al., 2014; Giannakeas et
al., 2011, 2012). When dealing with spin-polarized fermion systems, the
leading interaction term comes from a $p$-wave channel and this was confirmed
in cold atoms relatively early on by (Günter et al., 2005). One of the major
achievements with low-dimensional Fermi gases was the study of the spin-
imbalanced system by (Liao et al., 2010), effectively demonstrating that one
may use a cold-atom setup to simulate pairing phenomena in systems that have
an imbalance in spin population. This is related to the so-called FFLO pairing
phase and whether this may be stabilized by confining a system to a low-
dimensional geometry. For further discussion we refer to (Guan et al., 2013).

###  1.5 Review outline

The main part of our review starts with an introduction to the early models
and their exact solutions in Section 2. We then discuss single-component cold
atomic systems in Section 3 and multi-component systems in Section 4 with a
focus on collective modes, energy transfer, tunneling and quench dynamics in
small systems. Each section will present fundamental examples that are
targeted at new entrants to the field, and then provide a survey of more
recent results for researchers looking to update their knowledge.

In Section 5, we give an introductory overview of the physics of impurities
and polarons in low-dimensional systems under the assumption of strong
interactions. This provides both an invitation to the field for the newcomer,
as well as a survey of some of the recent results in both theory and
experiment. We then proceed to discuss the theory of spin-1 bosonic gases in
Section 6, giving both overviews of the state-of-the-art, and give recent
examples of mean-field studies to determine phase diagrams. We then proceed to
discuss control and dynamics in Section 7, where we review proposals on driven
dynamics and transport, as well as the emerging field of quantum control. The
latter is particularly intricate for few-body systems and is a topic of
extensive interest at present. We hope that the discussion here will give the
reader a starting point (and motivation) for further exploration.

A very recent topic of great interest is that of quantum droplets. These been
discussed in many studies in the past few years using a variety of different
methods, and the first experiments are now in place to gauge present
theoretical understanding. Here we provide an extensive review of recent
studies on droplet formation in Section 8 with an initial focus on static
properties in order to have a platform on which to discuss new proposals on
how to investigate their dynamics.

In the last two sections of our review, we discuss methods and the physics of
systems that display features due to interactions that are not short range.
Since the field of cold atomic few-body systems relies on a combination of
analytical insights and advanced numerical techniques, we provide in Section 9
an overview of some of the general techniques that are in use today. This
includes exact diagonalization, Bethe ansatz, DMRG, ML-MCTDHX and mappings of
strongly interacting systems to spin-chain descriptions. The goal is for the
reader to be able to obtain a quick overview of different options for
quantitative studies, as well as to provide the right references to more in-
depth technical discussions of each individual method and available resources
for their implementation.

After the discussion of the methods, we move on to systems that go beyond the
assumption of short-range interactions in Section 10 and discuss examples of
the physics that long-range interactions can add to few-body systems. This
includes p-wave interactions, particles with dipole moments and dipole-dipole
couplings, as well as combinations of atoms and ions in traps.

Finally, we end our review by a summary and give our view of some of the
exciting future perspective that this field holds, both in terms of the
exploration of fundamental physics at the border of few- and many-body systems
in low dimensions, as well as their potential use in quantum technologies in
the future.

##  2 Fundamental models and exact solutions

Long before one-dimensional atomic quantum systems could be realized in the
laboratory, theoretical models that would later become extremely relevant in
the field of ultracold atoms were already well developed in the context of
mathematical physics. This is due mainly to the fact that, given certain
conditions, interacting one-dimensional quantum systems can be exactly solved,
in the sense that the wave functions and the entire energy spectrum can be
obtained without approximations. One of the key ingredients for such an
accomplishment is a mathematical tool developed by Bethe in 1931 (Bethe,
1931), now called the coordinate Bethe ansatz: a trial wave function of a very
general character upon which a set of constraints (given by the particular
interacting model) is enforced.

Bethe’s original work employed this approach to solve the Heisenberg chain, a
quantum spin model of great importance in the field of condensed matter
physics and quantum magnetism. Later developments brought about a series of
works which made use of this same tool to exactly solve several other models.
The ones that would turn out to be most valuable for experiments with atoms at
low temperature are models of particles in the continuum with short-range
interactions that are modeled by a Dirac delta function, also known as contact
interactions. Among these models, the particular cases of the Lieb-Liniger
model for bosons and the Gaudin-Yang model for two-component fermions provide
examples of great impact that the applications of the Bethe ansatz would have
on current developments in physics (Batchelor, 2007). Comprehensive accounts
of both the bosonic and fermionic models described here, including details on
the Bethe ansatz solutions for each case, have been provided in (Cazalilla et
al., 2011; Guan et al., 2013).

We therefore discuss, here and in the following sections, the solutions and
general features of these models, with a focus on the few-body limit. In
Section 9, which is dedicated to a review of the underlying methods, we
provide a formal description of the solutions for arbitrary $N$. While we
focus on the coordinate version of the Bethe ansatz, a broader class of models
can also be exactly solved by other variants, such as the algebraic Bethe
ansatz (which is an essential ingredient of the celebrated Quantum Inverse
Scattering method (Takhtadzhan and Faddeev, 1979; Korepin et al., 1993;
Faddeev, 1995)), the analytical Bethe ansatz (Mezincescu and Nepomechie, 1992;
Foerster and Ragoucy, 2007) and the off-diagonal Bethe ansatz. (Wang et al.,
2015)

###  2.1 Bose gas with contact interactions

One of the most important applications of the coordinate Bethe ansatz for a
model in the continuum was put forward by Lieb and Liniger in 1963 (Lieb and
Liniger, 1963; Lieb, 1963). It addresses a bosonic ensemble of $N$ particles
with contact interactions in a one-dimensional system of length $L$ with
periodic boundary conditions. Nowadays, this model is especially relevant due
to its applications to the problem of interacting ultracold atoms. The
Hamiltonian for the Lieb-Liniger model reads

$H=-\sum_{i=1}^{N}\frac{\partial^{2}}{\partial
x_{i}^{2}}+2c\sum_{i<j}\delta(x_{i}-x_{j}),$ (2)

where $\hbar=m=1$ is assumed. The contact interactions between pairs of
particles are modelled by the delta-function term, with an interaction
strength determined by the parameter $c$, which can describe repulsive ($c>0$)
or attractive ($c<0$) regimes. It is customary to represent the interaction
strength in the Lieb-Liniger model in terms of the dimensionless parameter
$\gamma=c/\rho$, where $\rho=N/L$ is the particle density. Interestingly, in
contrast to the three-dimensional case, here we see that interactions grow
stronger as the density is reduced. More specifically, for fixed $c$ we find
that $\gamma$ becomes large at low density. This is an important feature that
needs to be kept in mind when working with one-dimensional interacting quantum
systems.

The Bethe ansatz wave function for this model is written as a linear
combination of $N!$ plane waves

$\psi(x_{1},...,x_{N})=\sum_{P}A(P)e^{i\sum_{j=1}^{N}k_{Pj}x_{j}},$ (3)

where $\\{k_{P1},...,k_{PN}\\}$ is a set of numbers called quasimomenta, or
sometimes rapidities. Here, $P$ describes a particular permutation of these
numbers and $A(P)$ denotes the permutation-dependent coefficients. In this
expression, it is assumed that the wave function describes a single coordinate
ordering, that is, $x_{1}<x_{2}<...<x_{N}$. The remaining sectors can be
obtained by simply considering the bosonic symmetry of the wave function upon
permutation of the particles positions. Enforcing the contact interaction
conditions upon the derivatives of the wave function, as well as the periodic
boundary conditions, results in a set of $N$ equations for $k_{j}$ which are
known as the Bethe ansatz equations and fully determine the solution of the
problem. In Section 9 we present details for these solutions.

###  2.2 Solutions in the few-body limit

To understand the nature and solutions of the Lieb-Liniger model for
interacting bosons, it is instructive to consider its simplest application,
namely the two-particle problem (Pethick and Smith, 2002). In this context,
Eq. (3) is given, in the region $0<x_{1}<x_{2}<L$ as
$\psi(x_{1},x_{2})=A_{12}e^{ik_{1}x_{1}}e^{ik_{2}x_{2}}+A_{21}e^{ik_{1}x_{2}}e^{ik_{2}x_{1}}$.
Assuming contact interactions leads to a relation for the coefficients
$A_{12}$ and $A_{21}$ of the form

$\frac{A_{21}}{A_{12}}=\frac{k_{2}-k_{1}+ic}{k_{2}-k_{1}-ic},$ (4)

where the bosonic symmetry upon permutation of particles has been enforced.
Minimizing the energy yields, for the ground state, $k_{1}=-k_{2}=k$, and the
addition of periodic boundary conditions allows us to write the resulting
ground state wave function for $N=2$ as

$\psi(x_{1},x_{2})=A\cos{\left[k\left(|x_{1}-x_{2}|-\frac{L}{2}\right)\right]},$
(5)

where $A$ is a normalization constant, and the quasimomentum $k$ and the
interaction strength $c$ are related by $c=2k\tan{\left(kL/2\right)}$. The
wave function in Eq. (5) respects bosonic symmetry when coordinates are
exchanged and vanishes for $x_{1}=x_{2}$ when $c\rightarrow\infty$, where we
also see that $k\to\pi/2$. Fig. 6 shows an illustration of the wave function
given in Eq. (5).

Figure 6: Ground state wave function of two interacting particles in a system
of length $L$ with periodic boundary conditions. Each curve belongs to a
different value of the interaction parameter $\gamma=c/\rho$. At
$\gamma\rightarrow\infty$, the wave function vanishes for $x_{1}=x_{2}$. Image
from (Pethick and Smith, 2002).

Formally, the three-particle bosonic problem is simply an extension of the
$N=2$ case described above. However, it is also an important example of a
minimal model where integrability is explicitly manifested (Jimbo, 1989). The
coefficients $A_{abc}$ corresponding to the wave function for $N=3$ respect
six relations analogous to that of Eq. (4) (Takahashi, 2005). They can be
transformed into one another through the action of
$Y_{ab}=-(k_{a}-k_{b}+ic)/(k_{a}-k_{b}-ic)$, a set of objects that in turn
obey

$Y_{ab}Y_{ac}Y_{bc}=Y_{bc}Y_{ac}Y_{ab}.$ (6)

Quite remarkably, this expression guarantees the existence of an exact
solution and also holds for the $N$-body problem. In the three-body regime, it
is the most elementary example of a Yang-Baxter equation (Yang and Ge, 1991;
Batchelor and Foerster, 2016), a relation that can be physically understood as
a pairwise factorization of interactions in the system. Since the operators
$Y$ permutate the quasimomentum values, the Yang-Baxter equation can be
interpreted as a condition of nondiffraction, meaning that the sets of
incoming and outgoing quasimomenta in scattering events are the same (McGuire,
1964; Sutherland, 2004). The Yang-Baxter relations are also respected by the
transfer matrices in two-dimensional lattice models (Baxter, 1972), and other
interpretations of this equation exist in different domains (see, for instance
(Vind et al., 2016; Batchelor and Foerster, 2016)). A set of operators obeying
the Yang-Baxter equation is generally a signature of integrability in many-
body systems (Batchelor, 2014). A precise definition of the concept of quantum
integrability, however, still remains a point of debate (Caux and Mossel,
2011), and one typically needs to treat systems on a case by case basis with
regards to questions of quantum integrability.

In the thermodynamic limit ($N,L\rightarrow\infty$ with the density $\rho=N/L$
fixed), the Bethe ansatz equations can be expressed as a Fredholm integral
equation for the momentum densities $f(k)$ (Takahashi, 2005). In Fig. 7, we
show Lieb and Liniger’s result for the energy of the one-dimensional bosonic
system with comparisons to known perturbative methods. Note that the
$\gamma\to\infty$ limit is approached very slowly in Fig. 7. This is the limit
where the bosonic system fermionizes as we will discuss in more detail below.
The thermodynamic properties of the Lieb-Liniger model at finite temperature
were obtained by Yang and Yang in 1969 with an approach that became known as
the thermodynamic Bethe ansatz (Yang and Yang, 1969), or the Yang-Yang
thermodynamic formalism. This method was later adapted to lattice problems,
such as the Heisenberg and Hubbard models (Takahashi, 1972, 2005; Essler et
al., 2005). Further developments of this approach include the Quantum Transfer
Matrix formalism (see, for instance (Suzuki, 1985; Klumper, 1992; Pâţu et al.,
2020)). While initially the main interest motivating the use of the Bethe
ansatz was finding the solutions and properties of many-body systems, its
validity extends to systems of any size, with arbitrary numbers of particles.
Indeed, solutions in the few-body limit can be found analytically and provide
great insight into the particular details of the method.

Figure 7: Comparison of the ground state energy of a repulsive Bose gas
obtained with a numerical solution of the Bethe ansatz equations to results
from perturbation theory. The dashed line at $e=\pi^{2}/3$ marks the result in
the limit of $\gamma\to\infty$. Image from (Lieb and Liniger, 1963)

###  2.3 The Tonks-Girardeau model

A particular limit of the Lieb-Liniger model of great importance for both
theory and experiment is the regime of infinitely repulsive interactions
$(\gamma\rightarrow\infty)$. This regime is characterised by surprisingly
simple solutions which are given by the so-called Bose-Fermi mapping, first
proposed by Girardeau in 1960 (Girardeau, 1960). It states that the wave
function $\psi_{B}$ for “impenetrable core" bosons can be written as

$\psi_{B}(x_{1},...,x_{N})=\prod_{i<j}\text{sgn}(x_{i}-x_{j})\psi_{F}(x_{1},...,x_{N}),$
(7)

where $\psi_{F}(x_{1},...,x_{N})$ is a wave function for identical spinless
fermions. The main feature of the mapping is to identify the bosonic wave
function as a symmetrized version of the fermionic wave function, which is
antisymmetric under permutation of particle coordinates. The latter is
obtained, for a system of $N$ particles, as the Slater determinant of single-
particle solutions $\phi_{n}(x_{i})$ for a given potential $V(x)$, where $n$
denotes the energy level. As a concrete example of this, the ground state
solution for the bosonic problem in a system of length $L$ with periodic
boundary conditions can be written, exploring the properties of the
determinant as222Formally, this solution holds only for odd $N$ in the
specific case of a ring geometry. For trapped systems, Girardeau’s mapping
holds for arbitrary $N$.

$\psi_{B}(x_{1},...,x_{N})=C\prod_{i<j}\Big{|}\sin{\left[\frac{\pi}{L}(x_{j}-x_{i})\right]}\Big{|}$
(8)

where $C$ is a normalization constant. The corresponding expression for a
harmonically trapped system has the form (Girardeau et al., 2001)

$\psi_{B}(x_{1},...,x_{N})=C^{\prime}\prod_{i=1}^{N}\exp\left(-\frac{x_{i}^{2}}{2b^{2}}\right)\prod_{i<j}|x_{j}-x_{i}|,$
(9)

where $b$ is the corresponding oscillator length, and $C^{\prime}$ is a
normalization constant. Note the products of coordinate differences in both
these cases. This structure is a result of the relation to Vandermonde
determinants of the antisymmetric expressions, and is for instance also found
in Laughlin wave functions at integer filling in the quantum Hall effect
(Scharf et al., 1994).

The Bose-Fermi mapping, which is restricted to one-dimensional geometries,
also guarantees that the energies and spatial distributions between the
bosonic wave function and its fermionic counterpart will be identical.
Furthermore, it provides a way of understanding how the strong interactions
between the bosons effectively mimic the Pauli exclusion principle, with the
wave function vanishing for $x_{i}=x_{j}$. The mapping does not hold, however,
for the momentum distribution, where the symmetry of the wave function plays a
significant role (Lenard, 1964).

In the attractive interaction regime ($c<0$), the ground state of $N$
particles is bound and unstable in the thermodynamic limit due to the scaling
of the energy of the system with the particle number. For finite systems, the
Bethe ansatz has been employed in the investigation of systems with periodic
(Sakmann et al., 2005a) and hard-wall (Hao et al., 2006b) boundary conditions,
again displaying bound configurations in the ground state which are
manifested, generally, as highly localized particle densities. Interestingly,
in the infinitely attractive limit ($c\rightarrow-\infty$), an excited state
with a gas-like configuration can be found. In analogy to the repulsive case,
it has been dubbed the Super Tonks-Girardeau regime, and it surprisingly
features properties which are similar to a strongly repulsive one-dimensional
gas (see Fig. 5), with even stronger correlations. This regime has been
accessed theoretically with Monte-Carlo simulations where the many-body wave
function is constructed from the exact solution of the two-body problem
(Astrakharchik et al., 2005). A Bethe ansatz treatment of the Super Tonks-
Girardeau gas has also been provided in (Batchelor et al., 2005).

###  2.4 Fermionic systems

The simplest example of an exactly solvable fermionic model is that of an
identical (or fully polarized) system of fermions. Due to the Pauli principle,
such a system is non-interacting and the Bethe ansatz need not be applied.
Instead, the wave functions for such systems are generally given in terms of
Slater determinants. In the multicomponent Fermi problem, however,
interactions between particles in different states are allowed. The Bethe
ansatz solutions to different versions of this problem have been provided
throughout the years in increasing degrees of complexity. Due to its general
character and applicability, a work that stands out in this context is Gaudin
and Yang’s solution to the case of a two-component fermionic system with an
arbitrary number of particles in each spin state
(${N_{\uparrow},N_{\downarrow}}$) (Gaudin, 1967; Yang, 1967). Two important
distinctions with respect to the previously explored bosonic case must be
highlighted: due to the Pauli principle, only interactions between fermions in
different spin states are allowed. Additionally, the resulting wave function
must respect the fermionic antisymmetry with respect to permutation of
identical particles. The solution for this Hamiltonian is given by a
generalized Bethe ansatz and a corresponding set of so-called nested Bethe
ansatz equations (see Section 9).

The simplest case of a fermionic system (excluding the spinless Fermi gas) is
the situation of a polarized system with a single opposite spin, that is
$N_{\downarrow}=1$. This problem was solved by McGuire both for repulsive
(McGuire, 1965) and attractive (McGuire, 1966) interactions. The situation
with two reversed spins, $N_{\downarrow}=2$ was further solved by Flicker and
Lieb (Flicker and Lieb, 1967). These studies show how reversing spins in a
fully polarised system can be energetically advantageous, in accordance with
the Lieb-Mattis theorem (Lieb and Mattis, 1962) stating that the ground state
should have the lowest possible total spin, i.e. for a balanced system the
total spin should vanish. The fermionic problem with higher spin was
investigated, with a focus on the repulsive regime, by Sutherland (Sutherland,
1968), while the attractive regime was later solved by Takahashi (Takahashi,
1970). The thermodynamic properties of the multicomponent fermionic gas were
obtained by Pedro Schlottmann (Schlottmann, 1993, 1994). The Bethe ansatz has
further been applied to the study of two-component bosons (Li et al., 2003;
Klümper and Pâţu, 2011; Klauser and Caux, 2011; Pâţu et al., 2018) and Bose-
Fermi mixtures (Lai and Yang, 1971; Imambekov and Demler, 2006b). The
thermodynamic properties of Bose-Fermi mixtures were obtained later with the
thermodynamic Bethe ansatz (Yin et al., 2009).

Two comments are in order regarding the applicability of the solutions for the
integrable systems described above: the first is the fact that, both at zero
and finite temperatures, static and dynamical properties of these models
assume the equilibrium condition. While extensions to non-equilibrium
situations (e.g. quantum quenches) have been explored in some studies (Zill et
al., 2016), this problem has remained a challenging one, inspiring development
of new approaches (see Section 9 for examples).

##  3 Single component static and dynamics

The ability to treat systems with small particle numbers exactly offers the
possibility to study the dynamics of not only classical properties, such as
density or momentum distributions, but also non-classical properties such as
coherence and entanglement. To be able to do this we will in the following
review different approaches for describing small systems, mostly in
inhomogeneous environments, and briefly discuss fundamental solutions and
equilibrium properties as a basis for understanding the dynamical properties.

###  3.1 Two-particles in a harmonic trap

The fundamental building block of many-particle systems is the two particle
system and a lot of ideas and concepts already emerge at this level. We
therefore begin with discussing the fundamental solution in 1D for two
particles of mass $m$ confined by a harmonic trapping potential of frequency
$\omega$ and interacting via a short-range delta-function interaction of
strength $g_{\text{1D}}$333Note that in free space this quantity relates to
the parameter $c$ in Eq. (2) as $g_{\text{1D}}=c\hbar^{2}/m$. (Busch et al.,
1998). In one dimension the Hamiltonian is given by

$H=\sum_{n=1}^{2}\left(-\frac{\hbar^{2}}{2m}\nabla_{n}^{2}+\frac{1}{2}m\omega^{2}x_{n}^{2}\right)+g_{\text{1D}}\delta(x_{1}-x_{2}),$
(10)

and due to the quadratic form of both the kinetic energy and the particle-
independent trapping potential it can be separated into a centre-of-mass and a
relative part, $H=H_{+}+H_{-}$, with

$\displaystyle H_{+}=$
$\displaystyle-\frac{1}{2}\nabla^{2}_{+}+\frac{1}{2}x_{+}^{2},$ (11)
$\displaystyle H_{-}=$
$\displaystyle-\frac{1}{2}\nabla^{2}_{-}+\frac{1}{2}x_{-}^{2}+g\delta(x_{-}).$
(12)

Here $x_{+}=\frac{x_{1}+x_{2}}{\sqrt{2}}$ is the centre-of-mass coordinate and
$x_{-}=\frac{x_{1}-x_{2}}{\sqrt{2}}$ is the relative coordinate. To achieve
the simple forms for $H_{+}$ and $H_{-}$, we have used harmonic oscillator
scaling where energies are given in units of $\hbar\omega$, lengths in units
of the ground state width $a=\sqrt{\hbar/m\omega}$ and the interaction
strength in units of $\sqrt{2\hbar^{3}\omega/m}$.

One can immediately see that the centre-of-mass Hamiltonian, $H_{+}$, simply
describes a single particle harmonic oscillator with the well known solutions

$\psi_{n}(x_{+})=(2^{n}n!)^{-1/2}\pi^{-1/4}\mathcal{H}_{n}(x_{+})e^{-\frac{x_{+}^{2}}{2}},$
(13)

where $\mathcal{H}_{n}$ is the $n$th order Hermite polynomial, and the
energies are given by $E_{n}=(n+1/2)$. One is therefore left with solving the
stationary Schrödinger equation for the relative Hamiltonian, which consists
of a single particle harmonic oscillator with a contact potential at the
origin. This immediately suggests that the odd eigenfunctions of the standard
harmonic oscillator are still eigenfunctions of $H_{-}$, as they possess a
node at $x_{-}=0$ and do not feel the additional contact potential. The even
states, however, are affected and can be found to be

$\phi_{n}(x_{-})=\mathcal{N}_{n}e^{-x_{-}^{2}/2}U\left(\frac{1}{4}-\frac{E_{n}}{2},\frac{1}{2},x_{-}^{2}\right),$
(14)

where the $U(a,b,z)$ are the Kummer functions and $\mathcal{N}_{n}$ is the
normalization. The corresponding eigenenergies are given by the implicit
equation

$-g=2\frac{\Gamma\left(-\frac{E_{n}}{2}+\frac{3}{4}\right)}{\Gamma\left(-\frac{E_{n}}{2}+\frac{1}{4}\right)}\;,$
(15)

and the spectrum is shown in Fig. 8. One can see that, as one would expect,
the odd eigenenergies are not affected by a change in $g$ and the even ones
increase and decrease depending on the sign of $g$. In particular, it is worth
noting that for $g\rightarrow\infty$ the energies of neighbouring odd and even
states become degenerate and the form of each even eigenfunction approaches
the absolute value of the next-higher-lying odd one (cf. Section 2.3).

Figure 8: Spectrum of the relative part, $H_{-}$, of the two-particle
Hamiltonian as a function of the interaction strength. The eigenvalues
corresponding to odd eigenstates are not affected (blue lines), whereas the
eigenvalue of the even eigenstates increase and decrease (red lines) for
larger positive and negative values of $g$.

While this model only contains essential elements, its applicability to
ultracold atom experiments was demonstrated by the Esslinger group (Stöferle
et al., 2006) and a large number of subsequent works have taken advantage of
it over the last two decades. More complex versions of the above model have
also been developed, for example considering anisotropic traps (Idziaszek and
Calarco, 2005; Bougas et al., 2020) or interactions via higher order partial
waves (Kanjilal and Blume, 2004; Idziaszek and Calarco, 2006; Zinner, 2012). A
further extension to the situation where the two particles are trapped in
differently centered harmonic oscillator potentials of the same frequency was
analytically solved by (Krych and Idziaszek, 2009), who showed that in this
case the relative Hamiltonian corresponds to a harmonic oscillator with an
off-center delta function. It is worth noting that many solutions to the
situations named above also exist in higher dimensions, but are not discussed
in this review.

It is also possible to extend the Hamiltonian in Eq. (10) to include the
effect of spin orbit coupling. The spectrum of such a system was discussed in
(Guan et al., 2014), the spin structure in (Guan and Blume, 2015) and a
numerical treatment that included the Raman coupling term was presented in
(Usui et al., 2020). Finally, the physics of two particles in an anharmonic
potential, where the center-of-mass and relative degrees of freedom can
couple, was explored in (Matthies et al., 2007), and for the situation of an
anharmonic optical lattice in (Mentink and Kokkelmans, 2009).

The existence of analytical solutions allows one to have full access to all
classical and quantum properties of the system, and we will discuss different
stationary properties below. However, it is worth noting that the existing
solutions can also be used to further applications in other areas. Examples of
this are a more precise description of the interaction term of the Bose-
Hubbard model (Büchler, 2010; Kremer et al., 2017), the description of quantum
gates (Calarco et al., 2000; Negretti et al., 2005), or for discussing the
measurement of spatial mode entanglement in ultracold bosons (Goold et al.,
2009). Furthermore, the knowledge of the exact solutions for the delta-split
harmonic oscillator has allowed one to study the many-body physics in the
Tonks limit in such systems (Busch and Huyet, 2003; Goold and Busch, 2008;
Goold et al., 2008).

Finally, the two-particle physics beyond the point-like interaction has been
discussed in the literature (see also Section 10.3). In (Ołdziejewski et al.,
2016) the physics of two particles with dipolar interaction was treated (in
full three-dimensions) and the stationary properties of two interacting
molecules were determined in (Dawid et al., 2018). Effects stemming from
combining anisotropic traps and realistic interaction potentials were
numerically discussed in (Grishkevich and Saenz, 2009).

###  3.2 Three particles

While systems with more than two particles in harmonic traps do not allow for
analytically exact solutions in all generality, considering the three particle
problem can still be very helpful. The Hamiltonian for the most general system
of two particles of species $A$ interacting with strength $g_{A}$, and one
particle of species $B$ which interacts with the $A$ particles through
$g_{AB}$, is given by

$H=\sum_{n=1}^{3}\left(-\frac{1}{2}\nabla_{n}^{2}+\frac{1}{2}x_{n}^{2}\right)+g_{A}\delta(x_{1}-x_{2})+g_{AB}\left[\delta(x_{1}-x_{3})+\delta(x_{2}-x_{3})\right],$
(16)

assuming all masses to be equal and using harmonic oscillator scaling. This
Hamiltonian is separable when moving to Jacobi coordinates

$\displaystyle X$ $\displaystyle=$ $\displaystyle(x_{1}-x_{2})/\sqrt{2}\,,$
(17) $\displaystyle Y$ $\displaystyle=$
$\displaystyle(x_{1}+x_{2})/\sqrt{6}-\sqrt{2/3}\,x_{3}\,,$ (18) $\displaystyle
Z$ $\displaystyle=$ $\displaystyle(x_{1}+x_{2}+x_{3})/\sqrt{3}\,,$ (19)

which allow one to write the Hamiltonian as separate centre-of-mass and
relative Hamiltonians, $H=H_{+}(Z)+H_{-}(X,Y)$. The center-of-mass part,
$H_{+}(Z)=-\frac{1}{2}\nabla_{Z}^{2}+\frac{1}{2}Z^{2}$, again just describes a
harmonically trapped single particle, while the relative part can be written
as

$\displaystyle H_{-}(X,Y)=$
$\displaystyle-\frac{1}{2}\left(\nabla_{X}^{2}+\nabla_{Y}^{2}\right)+\frac{1}{2}\left(X^{2}+Y^{2}\right)$
$\displaystyle+g_{A}\delta(X)+g_{AB}\left[\delta\left(-\frac{1}{2}X+\frac{\sqrt{3}}{2}Y\right)+\delta\left(-\frac{1}{2}X-\frac{\sqrt{3}}{2}Y\right)\right]\,.$
(20)

For $g_{A}=g_{AB}=g$ the eigenstates of $H_{-}(X,Y)$ possess $C_{6\nu}$
symmetry, allowing one to write the interaction part of the Hamiltonian in
cylindrical coordinates, $\rho=\sqrt{X^{2}+Y^{2}}$ and $\tan\phi=Y/X$, as

$H_{\text{int}}=\frac{g}{|\rho|}\sum_{j=1}^{6}\delta\left(\phi-\frac{2j-1}{6}\pi\right)\;.$
(21)

This symmetry can significantly simplify the calculation of eigenstates when
considering a truncated Hilbert space of non-interacting states or help in
creating relevant trial wave-functions for which the conditions from particle
exchange are simpler (Harshman, 2012; García-March et al., 2014). As a
benchmark for approximate treatments, highly accurate spectra for three
interacting bosons or fermions were determined in (D’Amico and Rontani, 2014).

Figure 9: Spectrum as a function of the interaction strengths for (a) three
indistinguishable bosons and (b) two indistinguishable bosons and a third
distinguishable particle. The axes are scaled in harmonic oscillator units.
Figure adapted from (García-March et al., 2016b).

In Fig. 9(a) the spectrum of three indistinguishable bosons is shown as a
function of the interaction strength, $g$. Compared to the two particle case
discussed above, the addition of an extra particle creates more degeneracies
in the spectrum at $g=0$ and $g\rightarrow\infty$. For the case of two
indistinguishable bosons of species $A$ and a third distinguishable particle
of species $B$ the spectrum is shown in Fig. 9(b). As there is no requirement
for symmetric scattering between the $A$ and the $B$ particles, the spectrum
for this situation contains more states compared to the panel (a), which
highlights the complexity that is introduced when multiple components are
considered.

###  3.3 Approximate approaches

Even though a general solution for one-dimensional gases with contact
interactions for any number of particles exist in free space (see Section
2.1), when going beyond three particles it is often helpful to consider
approximate techniques that allow to include external trapping potentials or
give easier access to specific quantities. In particular, for small systems it
is possible to consider many-body variational approaches and in the following
we will discuss a number of different ones.

In (Brouzos and Schmelcher, 2012a) this problem for bosons was approached
using the relative part of the Hamiltonian, after the center-of-mass part was
decoupled (see Fig. 10). All particle interactions are contained in the
relative Hamiltonian and the authors constructed a set of $N$-body analytic
wavefunctions, known as correlated pair wavefunctions (CPWF), from the exact
two-body solutions in terms of parabolic cylinder functions $D_{\mu}$ as

$\Psi_{\text{CP}}=C\prod_{i<j}^{P}D_{\mu}\left(\beta|x_{i}-x_{j}|\right),$
(22)

where $P=\frac{N(N-1)}{2}$ is the number of particle pairs and $C$ is a
normalization constant. Taking into account the boundary conditions from the
point-like interactions yields a transcendental equation for $\mu$

$\frac{g}{\beta}=-\frac{2^{3/2}\Gamma\left(\frac{1-\mu}{2}\right)}{\Gamma\left(-\frac{\mu}{2}\right)}\;,$
(23)

and while $\beta$ is known in the limiting cases of zero and infinite
strengths interactions, it can be treated variationally for finite $g$. This
approach leads to excellent agreement with exact numerical calculations for
intermediate interaction strengths (Barfknecht et al., 2015, 2016) and it was
later generalised to multi-component fermion systems (Brouzos and Schmelcher,
2013) and also shown to match quantum Monte-Carlo results for bosons in
different potential settings (Brouzos et al., 2013).

Figure 10: (a) Comparison of the total energies obtained numerically (marks)
and analytically (lines) via the CPWF as a function of the coupling strength
$g$. The numerical results for $N=3,4$ and $5$ are obtained by exact
diagonalization and for $N=6$ by the multi-configuration time-dependent
Hartree method. Panels (b)-(d) show the analytically and numerically obtained
one-body densities for different particle numbers and coupling strengths. All
quantities are in harmonic oscillator units. Figure adapted from (Brouzos and
Schmelcher, 2012a).

While most variational ansatz approaches rely on the Jastrow split of the
wavefunction into the non-interacting solution (center of mass part) and the
interacting part (relative part), they normally require the specific form of
the pair-correlation function to be known. However, in (Kościk et al., 2018) a
general pair correlation function of the form

$\varphi(x)=\frac{1}{\alpha}\left(1-\lambda e^{-\alpha|x|}\right)$ (24)

was suggested, which is independent of the exact form of the two-body
interaction. Here $\alpha=\sqrt{2/N}$ and $\lambda$ is chosen to fulfill the
discontinuity condition whenever two particles scatter. This is a significant
simplification that still gives good results while being numerically
inexpensive. A similar idea was used in (Kościk and Sowiński, 2020) to find
the many-body ground-state of interacting p-wave fermions in a one-dimensional
harmonic trap.

A different approach was presented in (Wilson et al., 2014), where a
wavefunction for a gas in a harmonic trap was geometrically constructed from
the known Bethe ansatz wavefunctions close to the trap centre where the
interactions dominate ($\lambda<\Lambda$) and a Gaussian decay at the trap
edges ($\lambda>\Lambda$)

$\Psi(\lambda,\vec{\theta})=\left\\{\begin{array}[]{cc}\psi_{B}\left(\vec{\kappa},\lambda,\vec{\theta}\right),&\lambda<\Lambda\\\
A(\vec{\theta})\exp\left(-\alpha(\vec{\theta})(\lambda^{2}-\Lambda^{2})\right),&\lambda>\Lambda\end{array}\right.,$
(25)

where $\lambda$ and $\vec{\theta}$ are the radial and angular components in
hyperspherical coordinates and both, the decay parameter
$\alpha(\vec{\theta})$ and the normalization $A(\vec{\theta})$, are fixed by
the boundary conditions at $\lambda=\Lambda$. To find the ground state energy
the position of the boundary $\Lambda$ and the quasi-momenta $k_{i}$ are used
as variational parameters. The authors showed that this approach works in the
attractive and the repulsive regime and can be highly accurate when compared
to numerical and exact benchmarks for two- and three-particle systems.

A variational approach for systems consisting of two indistinguishable
fermions and one distinguishable particle was presented in (Loft et al.,
2015). Since for the non-interacting limit ($g=0$) and for the impenetrable
one ($1/g\rightarrow 0$) one can find exact solutions for the factorized wave-
function in Jacobi hypercylindrical coordinates (Volosniev et al., 2015), the
authors assumed that this factorisation approximately holds for finite
interaction strengths as well. Based on this they constructed a variational
ansatz and showed that such an approach can describe the energies and the
spin-resolved densities of the system very well.

A straightforward approach that uses an interpolatory ansatz for many-body
systems in an inhomogeneous setting with finite strength interactions was
presented in (Andersen et al., 2016) and subsequently developed into a
systematic framework for both static and dynamic few-body studies (Lindgren et
al., 2020). Here the idea was to use a linear combination of known eigenstates
in the extreme limits of the interaction strength, $g\rightarrow 0$ and
$g\rightarrow\infty$, as the system at intermediate-strength interactions
contains a mixture of qualities from both limits. This technique was
subsequently also applied to a systems of fermions with equal and unequal
masses (Pȩcak et al., 2017), and for both cases good agreement to exact
results for up to six particles was shown. The unequal mass few-body states in
the limit $g\rightarrow\infty$ had been calculated around the same time using
semi-analytical techniques (Dehkharghani et al., 2016). This was later
generalized to larger systems of contact interacting systems with unequal mass
in harmonic traps, and used to infer the potential of both integrable, and
ergodic quantum chaotic behavior in such strongly interacting systems
(Harshman et al., 2017).

A model for few repulsively interacting particles trapped in a one-dimensional
harmonic well that included intrinsic $N$-body correlations was introduced in
(Andersen et al., 2017). The Hamiltonian was a hybrid of the Calogero and the
contact-interaction model and given by

$H_{=}\frac{1}{2}\sum_{i=1}^{N}\left(-\frac{\partial^{2}}{\partial
x_{i}^{2}}+x_{i}^{2}\right)+\frac{\sqrt{2}g}{\rho}\sum_{\langle
i,j\rangle}\delta(x_{i}-x_{j}),$ (26)

where $\rho$ is the relative hyper-radius for all $N$ particles. The two-body
interactions strength therefore depends on the position of all other particles
and the authors used this model to discuss insights into separability,
symmetry, and integrability.

For strongly interacting systems one can construct a mapping onto Heisenberg
spin chains (Volosniev et al., 2014a; Deuretzbacher et al., 2014; Levinsen et
al., 2015; Yang et al., 2015) (see also Section 9.5). The basic idea is that
for $g\rightarrow\infty$ the particles are arranged in a particular order, as
the wavefunction vanishes whenever two approach each other. However, for small
but finite values of $1/g$, neighboring particles can exchange their
positions, which effectively corresponds to nearest-neighbor spin
interactions. One can therefore describe the system with an effective
Hamiltonian that treats the particles as confined to a lattice for example as

$\mathcal{H}\simeq\frac{1}{g}\sum_{i=0}^{N-1}\left[J_{i}\mathbf{S}^{i}\mathbf{S}^{i+1}-\frac{1}{4}J_{i}\mathbf{I}\right]+E_{0}\mathbf{I}.$
(27)

Here $J_{i}$ is the nearest neighbour exchange constant, $\mathbf{S}^{i}$ the
spin operator at site $i$, and $\mathbf{I}$ is the unity operator. The rest of
the Hamiltonian is constructed such that the strongly interacting state has
the energy $E_{0}$. This approach allows one to calculate analytically
properties of general two-component systems in any external potential and can,
in principle, be extended to systems with $SU(N)$ symmetry. In Section 9.5 we
discuss this technique in detail.

###  3.4 Few-particle properties

Using the models described above and others existing in the literature, a
large amount of work has been carried out over the last two decades exploring
properties in many different regimes. While densities and momentum
distributions are often the first quantities to be analyzed, we focus in the
following on properties stemming from the additional degrees of freedom in
few-particle systems.

For example, probing the correlations in few-body systems is not an easy task
as it requires full reconstruction of the density matrix. However, recent
theoretical progress has shown that such a task is possible through the
measurement of position and momentum correlation functions when two atoms are
trapped in a double well potential (Bonneau et al., 2018). This approach was
successfully implemented for two interacting fermions in (Bergschneider et
al., 2019) and the density and momentum correlation functions found are shown
in Fig. 3. One can see that for attractive interactions double occupancy of
each well is observed, while for repulsive interactions it is suppressed with
each atom occupying different wells. In the non-interacting limit both
particles are delocalised over both wells. In the momentum distribution one
can see that the single particle coherence that is present in the non-
interacting case is suppressed in the presence of strong interactions, both
repulsive and attractive, while two-particle coherences are visible as
interference patterns along the antidiagonal and diagonal respectively. The
density matrix can be constructed from these momentum correlation functions
allowing for the measurement of entanglement entropy and lower bounds on the
concurrence. In (Becher et al., 2020) this technique was extended for systems
with up to three particles, allowing for measurement of identical particle
entanglement and the effect of antisymmetrization on correlations.

The properties of systems with up to six particles in a double well geometry
were explored in (Zöllner et al., 2006b) using a multiconfigurational approach
over the whole range of interaction strength, and the same authors also
subsequently studied the excitation spectrum (Zöllner et al., 2007).
Analytical results in the limit of two particles in a harmonic trap with a
tunable zero-ranged barrier at the trap center were presented in (Murphy et
al., 2007) and compared to a system with a finite width barrier in (Murphy and
McCann, 2008). The breakdown of a two-mode approximation for this situation
was explored in (Dobrzyniecki and Sowiński, 2016).

More detailed stationary states for a small number of harmonically trapped
bosons with attractive interaction were explored in (Tempfli et al., 2008),
and three distinct classes of states were found to emerge: (i) $N$-body bound
states, (ii) bound states of smaller fragments and (iii) gas-like states that
fermionize, that is, map to ideal fermions in the limit of infinite attraction
(cf. Section 2.3). For two-boson systems in oscillator potentials with various
powers of anharmonicity several measures for the entanglement were numerically
compared in (Okopińska and Kościk, 2009; Okopińska, 2010) and the non-local
properties of two interacting particles in separate traps at zero and at
finite temperatures were calculated in (Fogarty et al., 2011).

As briefly discussed in the introductory section, early landmark experiments
on the behaviour of small systems have been a strong motivator for the
theoretical works described here. These experiments include the fermionisation
of two distinguishable fermions in a 1D harmonic potential (Zürn et al.,
2012), where for diverging interaction strength the energy and square modulus
of the wave function for the two distinguishable particles becomes the same as
for a system of two noninteracting identical fermions (cf. Section 2.3). This
confirmed also the model of confinement-induced resonances predicted by
(Olshanii, 1998) and the exact results of (Busch et al., 1998) in a new
setting. The appearance of fermionisation was also confirmed to be compatible
with the full quasi-particle theory of tunneling developed in (Rontani, 2012).

###  3.5 Dynamical properties and quenches; time-dependent modulations

Because the evolution equations of small systems can be integrated exactly,
they are very interesting candidates to explore the non-equilibrium dynamics
following a quench from a well-defined initial state. With today’s existing
technologies, almost any part of the Hamiltonian can be quenched: kinetic
energies using artificial gauge fields, external potentials using
electromagnetic fields and interactions using Feshbach resonances.

One of the experimentally most common non-equilibrium situations is the
expansion of a quantum state into free space during a time-of-flight
measurement. In this situation any external trapping potential is
instantaneously switched off and the state is allowed to freely disperse. For
a one-dimensional Lieb-Liniger gas this situation was first discussed in
(Buljan et al., 2008), where the authors assumed that the gas was initially
localised in free space. Using exact solutions it was shown that the system
enters a strongly correlated regime during the free expansion in the 1D
direction, with the asymptotic form of the wave function having the
characteristic of a TG gas. This is a result of the drop in density during the
expansion, driving the system into the strongly interacting $\gamma\to\infty$
regime.

The problem of a Tonks–Girardeau gas released from a harmonic trap was solved
analytically in (Collura et al., 2013) and it was shown that in the
thermodynamic limit the reduced density matrix of any subsystem ultimately
converges to a generalized Gibbs ensemble as a consequence of the
integrability of the model. This integrability was also the reason behind one
of the most notable non-equilibrium experiments, namely the demonstration of a
quantum Newton cradle in a split Tonks-Girardeau gas (Kinoshita et al., 2006).
Here a gas of 40 to 250 87Rb atoms was observed to not noticeably thermalise
even after thousands of collisions (see Fig. 11). This experiment was later
carefully analysed by (van den Berg et al., 2016) who showed that the dynamics
happens on two time scales, with one corresponding to a rapid relaxation due
to the trap-insensitive contact interactions and a slower one governing the
in-trap periodic oscillations.

Figure 11: Density oscillations of a TG gas after in the initial state each
atom was in a superposition of $\pm 2\hbar k$ momentum. Image from (Kinoshita
et al., 2006)

One of the fundamental dynamical effects a quantum system can display is the
non-equilibrium tunneling between two states. For single particles or non-
interacting systems this is a textbook problem. However for interacting
particle system the presence of scattering leads to interesting new effects
and behaviours. An ab initio investigation based on the multi-configuration
time-dependent Hartree method (see Section 9.4) of the tunneling dynamics of a
few-atom system in a double well trap was presented by (Zöllner et al.,
2008b), covering the full range of interaction strengths. The authors found
that for weak interactions approximate Rabi oscillations can be observed,
whereas for larger interaction energies correlated pair tunneling emerges.
Close to the fermionization limit these fragment and lead to a rather involved
population dynamics in the two wells. More recently the behaviour of
interacting two and three particle systems in a quenched double well potential
was studied using a MCTDH-X method (Schäfer et al., 2020). Starting from the
ground state in a simple harmonic oscillator, the authors studied the
(quasi)-adiabatic vs. diabatic switching of a central barrier. They calculated
the evolution of the single-particle density matrix and the related von
Neumann entropy and found that diabatic switching leads to efficient energy
transfer between all particles across a large number of excited many-particle
states, while a quasi-adiabatic ramp only excites low lying states. Their
results were confirmed for systems of up to 10 particles. In the strongly
interacting limit, a spin model description of a few-body system with a
central barrier varying in strength with time was discussed by (Volosniev et
al., 2016), and a proposal for how to drive not just exchange of two particles
but also exchange of two pairs of particles via an engineered barrier dynamics
was presented.

As an alternative to the involved numerical solutions, (Dobrzyniecki and
Sowiński, 2018) treated the tunneling problem using a two-mode model that was
not just based on the lowest single-particle eigenstates of the external
potential, but used a basis of properly chosen effective wave functions. These
were found by decomposing the initial many-body state into the basis of exact
eigenstates of the many-body Hamiltonian, so that the interactions between
particles and the specifics of the initial many-body quantum state were both
taken into account. This approach was shown to significantly increase the
accuracy of the evolution for intermediate interactions.

Figure 12: Strengths of the dynamic structure function (DSF) for the one body
operator $\bar{F}=\sum_{i=1}^{N}x_{i}^{2}$ for a system of $2$ bosons with
different contact interaction strengths $g=0,1.75,10,30$. The DSF is shown for
the centre-of-mass (c.m.) component (red lines) which is independent of $g$,
while the relative components (blue and green lines) are distributed over many
excited states with its largest peak being the breathing mode. All quantities
are in harmonic oscillator (h.o.) units. Image from (Ledesma et al., 2019).

In an effort to compare the exact dynamics of two-particles to the standard
mean-field treatment of Bose-Einstein condensates, (Sowiński et al., 2010)
compared the evolution of the largest eigenvalue of the reduced single
particle density matrix (RSPDM) with the Gross-Pitaevskii equation (GPE)
dynamics and found that the difference between them signified the presence of
non-classical correlations in the system. A way to experimentally probe the
correlations in an interacting two-particle systems was presented in (Ledesma
et al., 2019), where the response to a monopolar excitation was studied via
the dynamical structure factor as a function of different values of the
interaction strength (see Fig. 12). This allowed the authors to compute
monopolar excitation energies, which are known to be directly related to the
internal structure of the system (Pethick and Smith, 2002).

Deeply connected to the question of coherence is the behaviour of two
interacting atoms in an interferometer setting and in (Fogarty et al., 2013)
it was shown that two interacting atoms impinging on a beam-splitter can
evolve into states close to the NOON-state (see Fig. 13. Creating entanglement
by scattering small numbers of atoms off each other in a harmonic traps was
investigated in (Holdaway et al., 2013). For four interacting particles the
authors were also able to identify parameter regimes in which attractive
interactions lead to a quantum soliton-like state, for which the time average
of the von Neumann entropy does not tend to a long-term mean value in contrast
to the case of repulsive interactions of similar magnitude. The idea to
control the motional entanglement between two colliding atoms by dynamically
varying the frequency of their common trap was explored in (Häberle and
Freyberger, 2014) and it was shown that the amount of entanglement can be
significantly increased with this dynamical tool.

Figure 13: Left: Single-particle density versus time for two attractively
interacting atoms. The barrier is point-like and positioned at $\tilde{x}=0$
and the two particles are initially trapped at $\tilde{x}=6$. The times for
scattering off the barrier ($t_{s}$), coming to rest at time at the classical
turning point ($t_{A}$), recombining and scattering a second time ($3t_{s}$),
and come to rest again ($t_{B}$) are indicated. Also shown is the Quantum
Fisher Information as a function of interaction strength $g$ and the barrier
height $\kappa$ at $t=t_{A}$ (middle panel) and $t=t_{B}$ (right panel). All
quantities are in harmonic oscillator units. Image adopted from (Fogarty et
al., 2013)

Since for the two-particle system in a harmonic trap the exact eigenstates are
known, it is interesting to see if analytical expressions can be obtained for
non-equilibrium quench dynamics as well. (Schmitz et al., 2013) explored the
breathing dynamics following a trap quench for a two-particle system and,
using the exact solution, calculated the full breathing mode spectrum. The
authors found that it consists of many frequencies which are accompanied by a
multiple of sidebands each. Using numerical methods for more particles, they
were also able to show how this case evolves towards the single-frequency
behavior that can be observed for larger particle numbers.

For a two-particle system (Budewig et al., 2019) studied an interaction quench
from repulsion to attraction and were able to derive analytical expressions
for the expansion coefficients of the time-evolved two-body wavefunction.
Quenching between zero to infinite interaction strength was discussed in
(Kehrberger et al., 2018), which is a situation in which analytical
expressions for the dynamical evolution of the wavefunction can be derived.
Quantifying the quench dynamics for small bosonic systems after a trap or
interaction change using thermodynamical quantities was studied in detail in
(García-March et al., 2016b) and a clear qualitative link between the amount
of irreversible work performed on the system and the increase in the degree of
inter-atomic entanglement was found. Similarly, in (Mikkelsen et al., 2022) it
was shown that the variance of the work probability distribution after a trap
quench is directly proportional to the infinite time-average of out of-time
order correlations functions (OTOCs) but only when finite interactions are
present between the particles.

A quantum quench from an ideal Bose condensate to the Lieb-Liniger model with
arbitrary attractive interactions was discussed in (Piroli et al., 2016a), and
in more detail in (Piroli et al., 2016b), where it was found that the
stationary state the system reaches at long evolution times after the quench
shows a prominent hierarchy of multi-particle bound states. While in
equilibrium the attractive Lieb-Liniger model is known to be thermodynamically
unstable, the appearance of these multi-particle bound states after the quench
ensures stability due to energy conservation. This is distinctively different
to the situation of the Super Tonks-Girardeau gas (cf. Section 2.3), which is
created by quenching from the strongly repulsive limit directly into the
strongly attractive regime. This state was first described by (Astrakharchik
et al., 2005) using a diffusion Monte Carlo method and corresponds to a stable
and highly excited Bethe state in the integrable Lieb-Liniger gas (Batchelor
et al., 2005; Chen et al., 2010). It was observed by (Haller et al., 2009),
who confirmed its stability due to the presence of strong quantum
correlations. As a step towards quantifying the stability of the Super Tonks
state, which is related to the rates of inelastic processes such as
photoassociation in pair collisions and three-body recombination, the local
two- and three-body correlators at zero and finite temperature were analyzed
in (Kormos et al., 2011).

While anyonic systems in the continuum have not yet been experimentally
realised, the equilibrium properties of the integrable anyonic Lieb-Liniger
model have been extensively studied (Pâţu et al., 2007; Batchelor et al.,
2007; del Campo, 2008), including also the case of small systems in harmonic
traps (Zinner, 2015). It was also shown how lattice setups can provide
realizations of the anyonic statistics known from the continuum models
(Bonkhoff et al., 2021). More recently, first steps towards studying the non-
equilibrium dynamics were taken. In (Piroli and Calabrese, 2017) the authors
considered a quench from the non-interacting ground state to the regime of
infinitely large repulsive interactions. While for anyonic systems one-point
functions and certain thermodynamic properties behave similarly to bosonic
systems (the anyonic parameter only renormalizes the pointwise interaction),
dramatic qualitative changes appear, for example, in the momentum
distribution: for a nonzero anyonic parameter it is asymmetric, indicating
that the Hamiltonian breaks parity symmetry. However, after the quench this
distribution evolves into a symmetric function for all value of the anyonic
parameter, and the system ‘forgets’ its anyonic character.

Studies focussing on periodically modulating the trapping frequency in 1D for
an interacting two- and few-particle systems were presented in (Brouzos and
Schmelcher, 2012b; Ebert et al., 2016) and a strong influence of the
interactions on the behaviour near parametric resonances was found. The effect
of a periodically oscillating magnetic field, which leads to a periodically
oscillating scattering length, was investigated for small fermion systems in
(Yin et al., 2016) and a technique to create controlled excitations based on
this effect was presented. The breaking of the scale invariance in a system
with smoothly varying external harmonic confinement in the nonadiabatic regime
through finite interparticle interactions was discussed in (Gharashi and
Blume, 2016).

Figure 14: Brody distribution parameter $\beta$ for $N=2,3,4$ interacting
bosons in a Kronig-Penney potential with $W$ wells as function of contact
interaction strength $\gamma$ and delta function barrier height $\tau$ between
adjacent wells. The Brody distribution is
$P_{\beta}(s)=(\beta+1)bs^{\beta}\exp(-bs^{\beta+1})$ with
$b=\left[\Gamma\left(\frac{\beta+2}{\beta+1}\right)\right]^{\beta+1}$ and
describes an integrable energy spectrum when $\beta\sim 0$ and a chaotic one
for $\beta\sim 1$. Black areas show $\beta<0$ indicating degeneracies and
picket-fence spectra. Image adopted from (Fogarty et al., 2021).

Similarly to quenching the interaction strength in a two-atom system, it is
interesting to consider quenching the presence of an impurity (see also
Section 5). In (Keller and Fogarty, 2016) the non-equilibrium dynamics
originating from such a process for an infinite mass impurity was considered,
using an approximate variational Lagrange-mesh method. Using the Loschmidt
echo to quantify the irreversibility of the system following the quench and
the probability distribution after long times to identify distinct dynamical
regimes, this work showed that the resulting dynamics depends strongly on the
interaction strength between the atoms. A similar situation, but allowing for
kinetic energy for the impurity, was considered in (Campbell et al., 2014),
where evidence of Anderson’s orthogonality catastrophe was found even for
small systems.

This situation was extended to a spinor impurity repulsively interacting with
a Bose-Einstein condensate after a quench in (Mistakidis et al., 2019c). Here
the temporal evolution was done using ML-MCTDHX and the determination of the
structure factor allowed the authors to identify three distinct dynamical
regimes as a function of increasing interspecies interactions. These regions
were found to be related to the segregated nature of the impurity and to the
Ohmic character of the bath. In the miscible regime one can identify the
polaron formation through the spectral response of the system, however for
increasing interaction an orthogonality catastrophe occurs and the polaron
picture breaks down. At this point the impurity transfers its energy to the
environment, which signals the presence of entanglement in the many-body
system. The situation for two interacting impurities was investigated in
(Mistakidis et al., 2020b), where the induced correlations between the
polarons due to the background were examined, see Section 5 for more details.

A further interesting question for quenched systems is the existence of
dynamical phase transitions. For the Tonks Girardeau gas this was investigated
in (Fogarty et al., 2017). Using the pinning transition for the situation
where a commensurate lattice was instantaneously turned on, this work
investigated the dynamics of the rate function and of local observables after
the quench. Whenever the return probability (or Loschmidt echo) goes to zero,
corresponding to singularities of the rate functions, the system is said to
display dynamical criticality and this work addressed the question of whether
this behaviour can be tied to the dynamics of physically relevant observables
or order parameters in the systems. Comparing the periodicity of the collapse
and revival cycle in the momentum distribution to that of the non-
analyticities in the rate function shows that the two are only closely related
for deep quenches.

Few-body systems have also been used to uncover the path from integrability to
chaos, as they allow one to explore equilibration and thermalization in
isolated systems. In (Fogarty et al., 2021) integrability breaking was
explored with interacting bosons trapped in multiwell potentials, and the
presence of both, finite interactions and finite barrier heights between the
wells, was shown to be able to result in strong signatures of quantum chaos,
which intensifies with increasing particle number. In Fig. 14 the integrable
and chaotic regions are discerned by fitting the energy spectrum to the Brody
distribution, with the fitting parameter $\beta\sim 0$ describing a Possonian
spectrum and therefore integrability, while $\beta\sim 1$ describes a Wigner-
Dyson distribution signifying strong energy level repulsion and therefore
chaos. Recent work by (Łydżba and Sowiński, 2021) has also investigated the
emergence of chaos between two component interacting fermions in a double well
potential and (Yampolsky et al., 2022) used a quantum generalization of the
Chirikov criterion of overlapping resonances to detect chaos in a two-particle
system in the presence of a perturbation. Finally, (Huber et al., 2021) used
convolutional neural networks to identify chaotic eigenstates of two fermions
and an impurity on a ring.

##  4 Two-component 1D gases

Multicomponent systems are not a straightforward extension of their single
component counterparts. Even small samples can already hold a large amount of
new physics, owing mostly to the presence of additional degrees of freedom,
which result in novel properties of the corresponding mixture. While single
component systems are defined by one scattering length, one external potential
and one dispersion relation, two component systems can have species-selective
trapping potentials, modified dispersion relations e.g. due to intercomponent
coupling (Pethick and Smith, 2002; Pitaevskii and Stringari, 2016), and two
additional scattering lengths, with one of them accounting for the inter-
species scattering. They therefore allow one to explore various correlation
regimes (Tempfli et al., 2009; García-March et al., 2014) and even feature an
experimentally realizable phase-separation phase transition (Timmermans, 1998;
Papp et al., 2008). Furthermore, one can also consider mixtures made from
particles of different statistics (Bose-Fermi mixtures), which we will address
below in Section 4.6.

Commonly, two-component bosonic settings consist of either two alkali metals
e.g. 23Na-87Rb, or two isotopes like 85Rb-87Rb, or two hyperfine states of the
same alkali metal (Myatt et al., 1997; Hall et al., 1998). The interplay
between the inter- and intraspecies scattering is a central issue in these
systems since it determines the static and dynamical properties of the mixture
and strikingly it allows for manipulating the degree of quantum correlations,
see Fig. 15. The role of the correlations has been extensively discussed in
the literature. Characteristic examples include, for instance, altered phase-
separation regimes (Garcia-March and Busch, 2013; García-March et al., 2014,
2014), the generation of quantum droplets (Petrov, 2015; Petrov and
Astrakharchik, 2016; Cabrera et al., 2018) (see in particular Section 8),
noticeable modifications of the superfluid to Mott insulator phase transition
(Łącki et al., 2013; Dutta et al., 2015), composite fermionization processes
(Zöllner et al., 2008a; Garcia-March et al., 2013; Zinner et al., 2014; Pyzh
et al., 2018), quantum emulsive states (Roscilde and Cirac, 2007), emergence
of instabilities (Cazalilla and Ho, 2003; Vidanović et al., 2013) and
intriguing phases such as paired or counterflow superfluidity (Hu et al.,
2009).

Figure 15: Schematic representation (the Miguel cube) of the few-body ground
states of a two-component bosonic mixture for different intra- and
intercomponent interactions. CF stands for composite fermionisation, PS refers
to phase-separation and FF stands for full fermionization. Figure from
(García-March et al., 2014).

###  4.1 Ground state properties of bosonic mixtures

The basic model for a pseudospinor 444Notice that multicomponent setting
containing $n$ individual components are often described as a fictitious
spin-$(n-1)/2$ system. They preserve the atom number of each component
prohibiting spin-mixing, while the internal states are not rotationally
symmetric in spin space, see also (Stamper-Kurn and Ueda, 2013). system of
ultracold atoms interacting via point-like interactions is described by the
Hamiltonian

$\displaystyle H$
$\displaystyle=\sum_{n=1}^{N_{A}}\left[-\frac{\hbar^{2}}{2m_{A}}\frac{\partial^{2}}{\partial
x_{n}^{2}}+V_{A}(x_{n})\right]+g_{A}\sum_{n=1}^{N_{A}}\sum_{m=n+1}^{N_{A}}\delta(x_{n}-x_{m})+\sum_{n=1}^{N_{B}}\left[-\frac{\hbar^{2}}{2m_{B}}\frac{\partial^{2}}{\partial
y_{n}^{2}}+V_{B}(y_{n})\right]$
$\displaystyle+g_{B}\sum_{n=1}^{N_{B}}\sum_{m=n+1}^{N_{B}}\delta(y_{n}-y_{m})+g_{AB}\sum_{n=1}^{N_{A}}\sum_{m=1}^{N_{B}}\delta(x_{n}-y_{m}).$
(28)

Here, the components $A$ and $B$ are under the influence of the external
trapping potentials $V_{A}(x)$ and $V_{B}(y)$ respectively. Such a setup can
be experimentally realized either by considering two distinct hyperfine states
of the same element (mass-balanced system, $m_{A}=m_{B}\equiv m$) or two
different atomic isotopes (mass-imbalanced setting, $m_{A}\neq m_{B}$). In the
former case the pseudospin degrees-of-freedom can correspond, for instance, to
the hyperfine states $\ket{A}\equiv\ket{F=1,m_{F}=1}$ and
$\ket{B}\equiv\ket{F=1,m_{F}=-1}$ of a 87Rb gas and the latter scenario
describes a mixture of isotopes containing e.g. 87Rb and 7Li. The effective 1D
intra-(intercomponent) coupling constants $g_{\sigma}$
($g_{\sigma\sigma^{\prime}}$) with $\sigma^{\prime}\neq\sigma=A,B$ are related
to the 3D $s$-wave scattering lengths, ${a^{s}_{\sigma\sigma^{\prime}}}$
(Olshanii, 1998). Experimentally, $g_{\sigma\sigma^{\prime}}$ can be tuned
either with the aid of ${a^{s}_{\sigma\sigma^{\prime}}}$ utilizing Feshbach
resonances (Köhler et al., 2006; Chin et al., 2010) or by means of the
transversal confinement frequency ${{\omega_{\bot}}}$ using the so-called
confinement-induced resonances (Olshanii, 1998), see also Section 10.1. Recent
advances in optical trapping techniques allow one to imprint almost arbitrary
shapes of the external trap and even create species-dependent potentials e.g.
by employing the so-called “tune-out” wavelength (LeBlanc and Thywissen, 2007;
Henderson et al., 2009).

Recently this topic has been comprehensively summarised in a review by
Sowinski and García-March (Sowiński and García-March, 2019), which mostly
focused on the stationary properties of few-body mixtures in harmonic traps. A
multitude of distinct phases can be formed stemming from the interplay between
intra- and intercomponent couplings, see Fig. 15 for an overview. Indeed,
there are six individual configurations for repulsive couplings (García-March
et al., 2014).

The so-called BEC-BEC limit appears when the particles are non-interacting and
the BEC-TG phase occurs for infinite intracomponent coupling of one species
and all other interactions vanishing. Systems enter the TG-TG phase if both
intracomponent interactions are infinitely strong and the intercomponent
coupling is zero. In the composite fermionization regime the intracomponent
couplings are suppressed, but strong intercomponent interaction exist.
Finally, the so-called phase-separation regime, is realized for infinite
intercomponent repulsion, while one of the intracomponent strengths has a
finite value and the other is zero. This state is also generated for finite
interactions as long as the intercomponent coupling is larger than the mean of
the intracomponent ones. In the finite interaction regime and in particular
for two different hyperfine states of 87Rb phase-separation and associated
mixing-demixing processes have been experimentally realized in quasi-1D traps
(Tojo et al., 2010; Eto et al., 2016a, b).

It is worth mentioning that phase separation occurs under the same conditions
in lattice two-component bosonic mixtures (Mishra et al., 2007; Alon et al.,
2006; Richaud et al., 2019) and its length scale (determined in terms of a
single lattice site) depends crucially on the value of the intraspecies
repulsions. The later determine the superfluid or Mott-insulator character of
the mixture. The underlying phase transition boundary is shifted when compared
to the single component scenario due to intercomponent entanglement as was
discussed for higher dimensions e.g. in (Wang et al., 2016; Buonsante et al.,
2008; Gadway et al., 2010). For a more elaborate discussion on the various
emerging phases and admixtures of energetically higher-bands via extended
Bose-Hubbard models, which are beyond the scope of the present review,
describing these systems see (Altman et al., 2003; Isacsson et al., 2005).

Coming back to Fig. 15, the full fermionization regime appears when all
interaction strengths tend to infinity. In the case of intercomponent
attractions and strong intracomponent repulsions it is possible to form a
molecular TG state among the species if the attraction is smaller or of the
same order as the repulsion (Tempfli et al., 2009). A further increase of the
attraction leads to a prominent bound state and eventually to the collapse of
the gas. This is discussed in further detail in Section 8. Finally, for strong
attractive intracomponent and vanishing intercomponent couplings the so-called
Super Tonks (excited) state occurs (see also Section 2.3), while tuning the
intercomponent attraction to larger values a corresponding molecular Super
Tonks phase can appear.

There are several ways to induce nonequilibrium quantum dynamics of
pseudospinor systems including, for instance, sudden changes or time-dependent
variations of the interatomic interactions and the external confinement.
Accordingly, it is possible to dynamically trigger complex tunneling
mechanisms beyond the single component paradigm, create peculiar correlation
patterns, realize few-body analogs of macroscopic phenomena, and trigger
collective or local excitations of the binary system to name a few. Below, we
review the observations of such processes in few-body and mesoscopic settings.

Analysis of pseudospinor systems requires theoretical methods and numerical
approaches that are able to account for both inter- and intracomponent
correlations, see Section 9. The latter are naturally enhanced for few-body
systems especially in 1D. To have a comparison point, let us first mention
that in the limit of large particle numbers and weak interactions, where
quantum fluctuations can be neglected, the standard approach is the celebrated
mean-field approximation. It is based on the assumption that each of the
species, $A$ and $B$, occupies only one single-particle state and can be
described by a macroscopic wave function $\phi_{A}(x)$ and $\phi_{B}(x)$,
respectively (Pethick and Smith, 2002; Kevrekidis et al., 2007). This
assumption reduces the many-body problem to an effective single-particle one,
where all interaction effects solely manifest themselves in the deformation of
the associated single-particle orbital.

The mean-field ansatz for the ground-state wave function according to the
variational principle leads to the coupled Gross-Pitaevskii equations

$i\hbar\frac{\partial\Psi_{A}(x,t)}{\partial
t}=\bigg{[}-\frac{\hbar^{2}}{2m_{A}}\frac{\partial^{2}}{\partial
x^{2}}+V_{A}(x)+g_{A}\left|\Psi_{A}(x,t)\right|^{2}+g_{AB}\left|\Psi_{B}(x,t)^{2}\right|\bigg{]}\Psi_{A}(x,t),$
(29) $i\hbar\frac{\partial\Psi_{B}(x,t)}{\partial
t}=\bigg{[}-\frac{\hbar^{2}}{2m_{B}}\frac{\partial^{2}}{\partial
x^{2}}+V_{B}(x)+g_{B}\left|\Psi_{B}(x,t)\right|^{2}+g_{AB}\left|\Psi_{A}(x,t)^{2}\right|\bigg{]}\Psi_{B}(x,t).$
(30)

Here each component is also subjected to an external potential $V_{i}(x)$ with
$i=A,B$ and it is assumed that $\int dx\left|\Psi_{i}(x)\right|^{2}=1$. The
intra and interspecies coupling strengths are denoted by $g_{i}$ and $g_{AB}$
respectively with $a_{i}$ and $a_{AB}$ being the corresponding scattering
lengths [see also Eq. (1)]. The Gross-Pitaevskii framework has been proven to
accurately describe the physics of weakly interacting single- and
multicomponent bosonic gases. It can explain a plethora of stationary and
nonequilibrium phenomena, e.g. the formation, interactions and dynamics of
various types of nonlinear excitations (Kevrekidis et al., 2007, 2015;
Kevrekidis and Frantzeskakis, 2016). In this review, such processes will not
be described in detail, as they usually do not have a well-defined few-body
counterpart. The focus is rather placed on processes where correlation-induced
phenomena play a dominant role.

In the so-called Thomas-Fermi limit where the interaction energy dominates,
the kinetic energy of the mixture can be neglected (Pitaevskii and Stringari,
2016; Kevrekidis et al., 2007). The condition for phase-separation then only
depends on the relative interaction strengths between the atoms and in
homogeneous systems the separation transition happens when
$g_{AB}>\sqrt{g_{A}g_{B}}$ (Ho and Shenoy, 1996; Ao and Chui, 1998; Pu and
Bigelow, 1998; Timmermans, 1998). For condensates in inhomogeneous potentials
usually regions exist where the kinetic energy of the atoms cannot be
neglected and, as a consequence parameters such as the trap frequency, atom
number and mass ratio can noticeably impact the miscible-immiscible phase
boundary (Wen et al., 2020; Tanatar and Erkan, 2000; Lee et al., 2016;
Gutierrez et al., 2021; Pyzh and Schmelcher, 2020). As an example, for mass-
imbalanced mixtures the degree of spatial localization of each component is
naturally affected by the unavoidable change of the related confinement length
scale (Pflanzer et al., 2009, 2010; Cikojević et al., 2018), see also (Pęcak
and Sowiński, 2019; Pȩcak et al., 2016; Cui and Ho, 2013; Harshman et al.,
2017; Volosniev, 2017) for related effects in 1D fermionic mixtures. In the
presence of a harmonic trap, the heavier component is usually stronger
localised and has a ground state density that has a Gaussian distribution of
width $\sqrt{\mathcal{\kappa}/2}$, where $\mathcal{\kappa}$ is the mass ratio.
It approaches a $\delta$-like shape for $\mathcal{\kappa}\to 0$. It is also
noteworthy that keeping the interspecies repulsion fixed and increasing solely
the mass ratio, one can also achieve a transition from a miscible to an
immiscible phase. This transition behavior can be seen on the one-body density
level and it is accompanied by a minimum in the two-body interspecies
correlation function, implying the formation of anti-ferromagnetic order
(García-March et al., 2014; Garcia-March and Busch, 2013). Moreover, for
strongly mass-imbalanced harmonically trapped few-boson mixtures and in the
limit of large intercomponent repulsions it has been shown in (Pflanzer et
al., 2009, 2010) that the heavier component acts as an effective “material
barrier” for the lighter species, inducing an effective attraction between the
atoms of the latter that increases with decreasing mass ratio. This behavior
can be explained by constructing an effective Bose-Hubbard model for the light
atoms which contains a renormalized hopping and on-site interaction
parameters.

###  4.2 Energy transfer and breathing modes

A fundamental mechanism that arises in multicomponent settings due to the
intercomponent coupling is dynamical "correlation transfer" between the
components. A detailed study of such process was presented in (Krönke et al.,
2015) with the focus on quantum dynamics of a single atom of species $A$ in a
finite bosonic environment made of $B$ atoms. Initially, both components are
trapped in species-selective 1D harmonic traps such that their spatial overlap
is rather small. Upon changing the initial displacement of the harmonic traps,
the two species collide (see e.g. (Vogels et al., 2002; Perrin et al., 2007)
for relevant experiments) and thereby interchange energy and correlations as
well as build up interspecies correlations (Deuar and Drummond, 2007; Ziń and
Wasak, 2018; Ögren and Kheruntsyan, 2009). In particular, the atom $A$
oscillates in the harmonic trap and thereby periodically penetrates the $B$
species held in a displaced trap. It has been shown that the corresponding
interspecies energy transfer accelerates with increasing the number of $B$
atoms which was related to a level splitting of the involved excited many-body
states. The finite amount of interspecies correlations created stimulates
coherent energy transfer processes which manifest in singlet and doublet
excitations of the bosonic environment.

Figure 16: Time-evolution of the density (in the relative frames) depicting
the breathing mode dynamics of a $N_{A}=N_{B}=2$ bosonic mixture following a
change of the trap frequency. $R_{AB}$ denotes the relative coordinate between
the individual center-of-masses of each species, while $r_{A}$ and $r_{B}$ are
the relative coordinates of the $A$ and $B$ species, respectively. Panels
(a)-(c) [(d)-(f)] represent the decoupled $g_{AB}=0$ [strongly coupled,
$g_{AB}=10$] interaction regime. The ensuing oscillation frequencies
correspond to energy differences between specific eigenstates of the system as
explained in (Pyzh et al., 2018). All quantities are expressed in postquench
harmonic oscillator units. Figure from (Pyzh et al., 2018).

The correlated breathing mode dynamics of few-body binary bosonic mixtures
confined in a 1D harmonic trap was investigated in (Pyzh et al., 2018). Using
exact diagonalization with a correlated basis approach, the low-energy
eigenspectrum was obtained for the complete range of repulsive intra- and
intercomponent interaction strengths. In particular, the strongly coupled
limiting cases of composite fermionization, full fermionization and phase-
separation as well as all the intermediate symmetric and asymmetric
intraspecies interaction regimes were examined in detail. Knowing the
stationary spectrum of the system, one can see that the breathing dynamics in
the linear-response regime can be triggered by slightly quenching the trap
frequency symmetrically for both components. Depending on the choice of
interaction strengths, the presence of one to three monopole modes was
identified besides the breathing mode of the center-of-mass coordinate, see
Fig. 16. In particular, it was shown that for the uncoupled mixture each
monopole mode corresponds to the breathing oscillation of a specific relative
coordinate. Importantly, increasing the intercomponent coupling leads to
multi-mode oscillations in each relative coordinate, which turn into single-
mode oscillations of the same frequency in the composite-fermionization
regime.

###  4.3 Quench dynamics and pattern formation

The interspecies interaction quench-induced phase-separation dynamics of a
harmonically trapped and repulsively interacting binary bosonic ensemble was
discussed in (Mistakidis et al., 2018). It was shown that the increase of the
interspecies repulsion results in a filamentation of the density of each
species. The spontaneously generated filaments are strongly correlated and
exhibit domain-wall structures at the two-body level as shown in Fig. 17.
Following the reverse quench protocol, i.e. decreasing the interspecies
interaction strength, the formation of multiple dark-antidark solitary waves
was observed (Katsimiga et al., 2020; Kiehn et al., 2019). The experimental
realization of such soliton configurations in a 87Rb BEC containing
$~{}10^{6}$ atoms was also recently reported (Katsimiga et al., 2020). These
solitary wave structures have been found to decay into the many-body
environment soon after their generation, which is in sharp contrast to the
predictions of the mean-field approximation. It constitutes a basic
characteristic of many-body soliton entities whose properties have been
extensively discussed, e.g. in (Katsimiga et al., 2017a, 2018; Dziarmaga and
Sacha, 2002; Mishmash and Carr, 2009; Delande and Sacha, 2014; Syrwid and
Sacha, 2017; Syrwid, 2021). Moreover, experimental counterparts of these
observations were reported by simulating single-shot images showcasing, among
others, that the growth rate of the corresponding variance of a sample of
single-shots555Detailed discussions regarding the theoretical implementation
of the single-shot procedure in MCTDH type methods can be found in (Sakmann
and Kasevich, 2016; Lode and Bruder, 2017; Katsimiga et al., 2017b, 2018) for
single component atomic ensembles and in (Mistakidis et al., 2018, 2019a;
Erdmann et al., 2019, 2018; Koutentakis et al., 2019) for multicomponent ones.
Note also that in single component settings the behavior of the single-shot
variance is related to the degree of fragmentation in the system. probes
entanglement in the system (Mistakidis et al., 2018; Katsimiga et al., 2017a,
b).

Figure 17: (Upper panels) Density evolution of the $(a)$ A and $(b)$ B species
following an interaction quench of a decoupled binary mixture into the
immiscible phase. (Lower panels) Snapshots of the two-body intracomponent
correlation function of the $(c)$ [$(d)$] A [B] species and $(e)$ the
interspecies two-body correlation function representing the domain-wall
formation among the species. Figures adapted from (Mistakidis et al., 2018).

The defect scaling and its relation to the celebrated Kibble-Zurek
mechanism666The Kibble-Zurek mechanism (Kibble, 1976; Zurek, 1985) is based on
the assumption that every system undergoing a phase transition experiences a
phase of “impulsive” behavior where the topological defects “crystallize” in
the system. It predicts the scaling of the density of the generated
topological defects with the rate at which the phase transition critical point
is crossed. when dynamically crossing the phase-separation boundary of binary
Rabi-coupled bosonic mixtures was studied in (Sabbatini et al., 2011, 2012)
within a truncated Wigner approximation. The latter method accounts for
quantum fluctuations around the mean-field solution (Blakie et al., 2008).
Adding a Rabi coupling term with time-dependent strength $\Omega(t)$ to Eq.
(28) will allow particle transfer among the components. Importantly, the
miscible to immiscible phase transition can be realized by appropriately
varying the Rabi amplitude above/below a critical strength $\Omega_{cr}$ in
addition to choosing a combination of the interaction values. A detailed
characterization of the ground state magnetic properties and excitation
spectrum of Rabi-coupled two-component bosonic systems 777The presence of
Rabi-coupling triggers spin-mixing processes among the internal states and
favors their miscibility Recati and Stringari (2021). Exemplary, experimental
demonstrations of such gases can be found in Williams et al. (2000); Zibold et
al. (2010); Nicklas et al. (2015); Cominotti et al. (2021). can be found in
(Abad and Recati, 2013; Sartori and Recati, 2013; Tommasini et al., 2003;
Recati and Stringari, 2021). Specifically, in (Sabbatini et al., 2011, 2012),
a mixture of $N\approx 10^{5}$ atoms was initially prepared in the immiscible
phase with $\Omega(0)>\Omega_{cr}$. Subsequently, $\Omega(t)$ was linearly
ramped-down to zero with a finite ramp rate $\tau_{Q}$ giving rise to long-
lived domain-wall configurations among the components. Remarkably it was shown
that the number of spontaneously generated domains in a ring trap of length
$L$ scales for varying quench rate with an exponent as predicted by the
Kibble-Zurek mechanism. In particular the number of defects according to
Kibble-Zurek theory is given by

$N_{d}=\frac{L}{\xi_{0}}\left(\frac{\tau_{0}}{\tau_{Q}}\right)^{1/3},$ (31)

where $\tau_{0}=\hbar/(2g_{S}\rho)$, $g_{S}=(g_{A}+g_{B}-2g_{AB})/2$ and
$\xi_{0}=\xi_{S}/\sqrt{2}$ is related to the spin healing length
$\xi_{S}=\hbar/\sqrt{2m\rho g_{S}}$, while $\rho$ denotes the uniform density
of the gas. In the actual numerical simulations it was observed that for
$\tau_{Q}>2$ms it holds that $N_{d}\propto\tau_{Q}^{-n}$ with $n=0.341\pm
0.006$ which is in good agreement with the exponent in Eq. (31). Further
deviations from this scaling law are present in elongated harmonic traps due
to the spatial inhomogeneity induced by the external confinement.

Figure 18: (a) Experimental setup for a 1D double-well potential created by
focusing two laser beams with a high-resolution objective. An acousto-optic
deflector (AOD) is utilized to independently control the intensity and
position of the two laser beams, which allows one to adjust the tunnel
coupling $J$ and the tilt $\Delta$ between the two wells. (b) Single-particle
(green line) and pair tunneling (blue line) probabilities at the interaction
energy $U/J=10.05\pm 0.19$ with respect to the tilt $\Delta$. The experimental
data points have $1\sigma$ statistical uncertainty (error bars) and represent
the time-averaged probability of detecting a single or a pair of particles in
the right well. The system is initially prepared with two atoms in the left
well and the the tunneling is switched on. As shown, pair tunneling is
resonant in a symmetric double-well, while conditional single-particle
tunneling takes place for $-2\Delta=U$. Figures adapted from (Murmann et al.,
2015a)

###  4.4 Tunneling dynamics

As discussed in Section 3 the tunneling dynamics of single species bosons in
the crossover from weak to strong interactions shows intriguing phenomena
including common Josephson oscillations, pair tunneling, self-trapping, and
fermionization (Zöllner et al., 2008b). Extensions of these studies to binary
bosonic mixtures were introduced in several works, e.g., in (Mazzarella et
al., 2011, 2009; Gillet et al., 2014), demonstrating additional complexity due
to the involvement of intra- and intercomponent interactions, as well as the
mass-imbalance between the components. Indeed, in such non-linear systems time
evolution cannot be simply explained by the knowledge of the single component
system except for the two cases of either $g_{AB}=0$ and $g_{A}=g_{B}$ or
$g_{AB}\to g_{A}$ with $g_{B}=0$ and $N_{B}=1$ in mass-balanced settings.
Knowledge of the structure of the full few-body spectrum is required.

Relevant dynamical aspects observed also for macroscopic bosonic ensembles
include, for instance, macroscopic quantum self-trapping and coherent quantum
tunneling (Kuang and Ouyang, 2000; Smerzi et al., 1997; Raghavan et al.,
1999), collapse and revival of population dynamics (Naddeo and Citro, 2010) as
well as symmetry breaking and restoring scenarios (Satija et al., 2009). The
Josephson or coherent quantum tunneling is a direct generalization of its
single component counterpart (see also Section 3.5), where the condensates
oscillate symmetrically among the two potential wells, thus being reminiscent
of the sinusoidal Josephson effect appearing in superconductors. Moreover, the
macroscopic quantum self-trapping refers to an essentially “locked” population
imbalance among the components. It is a quantum effect maintained in closed
systems that originates from the presence of non-linear interactions. However,
it strongly depends on the system parameters such as the initial state and the
number of atoms. Particularly, depending on the initial condition there are
two different symmetry-broken phases of macroscopic quantum self-trapping in
which the components reside either in distinct wells or at the same site.
Moreover, symmetry restoring tunneling modes exist, where the two species
either avoid or follow each other in the course of the evolution. Along the
same lines collapse and revival events which occur in the dynamics of the
population imbalance between the two species are pure interaction
effects888Precisely, both inter- and intracomponent interactions produce
collapses and revivals of the population imbalance which counteract each other
and are suppressed in the equal interaction limit. and constitute a
macroscopic quantum phenomenon. These phenomena are well-known in non-linear
optics and are captured by the Jaynes-Cummings model (Scully and Zubairy,
1999), which describes the interaction of a single-mode quantized radiation
field and a two-level atom.

The above tunneling processes have been observed in the large particle number
limit and for weak interactions where the relevant dynamics can be described
using the coupled set of mean-field Gross-Pitaevskii equations of motion given
in Eqs. (29) and (30) Raghavan et al. (1999); Satija et al. (2009), or using
the lowest-band approximation and especially the Bose-Hubbard model Kuang and
Ouyang (2000); Smerzi et al. (1997); Naddeo and Citro (2010). Note that the
particular case of the two-site Bose Hubbard model, which is one of the
simplest and most studied models that accounts for tunneling and self-trapping
in double wells, can be solved exactly by Bethe ansatz methods (Links et al.,
2006; Tonel et al., 2005). However, for few-body systems and strong
interactions, correlation phenomena become important and invalidate the
aforementioned models. Indeed, the study of few-body systems is particularly
useful, not only due to their numerical feasibility for ab-initio
calculations, but most importantly since they allow for an in depth
understanding of the accompanying microscopic mechanisms which subsequently
provide a bottom-up perspective to the processes occurring in larger systems.
Additionally, recent experimental developments enable one to design and probe
the tunneling properties of few-body systems. A characteristic example is the
experimental realization of the tunneling of two distinguishable fermions
where the interaction strength and the shape of the external potential were
used to engineer the process (Murmann et al., 2015a), see Fig. 18.

The role of intra- and intercomponent interactions in the tunneling dynamics
of few-body bosonic mixtures in a 1D double-well was discussed in (Chatterjee
et al., 2012) for different initial configurations, namely a population
imbalance and a phase-separated state. To elucidate the underlying few-body
correlation processes a mixure of three $A$-species bosons and a single
$B$-species atom was considered. Relying on the energy spectrum of the system,
it was found that stronger interspecies repulsions result in a larger
tunneling period. This mechanism is the few-body equivalent of quantum self-
trapping (Smerzi et al., 1997; Raghavan et al., 1999; Juliá-Díaz et al., 2010;
Wüster et al., 2012) which has been experimentally observed by (Albiez et al.,
2005) using 87Rb condensates, see also the review by (Gati and Oberthaler,
2007) and references therein. Therefore, the interspecies repulsion leads to
particle “binding” and reduces the tunneling rate. Moreover, the tunneling
behavior crucially depends on the initial configuration. Indeed, for complete
population imbalance corresponding to either all atoms prepared in the same
well or experiencing intercomponent separation, the individual components
tunnel in- or out-of-phase. However, for “partial” population imbalance where
one component is delocalized equally over the two wells and the other is
localized a counterflow type dynamics takes place. We note that these
oscillations among wells can also be observed in the limit of very strong
interactions (Barfknecht et al., 2018b).

Figure 19: Characteristic tunneling regions of $N_{A}=2$ bosons through a
single $B$ atom as captured by the population of the right well, see text for
details. Shown are different mass ratios, namely $m_{A}/m_{B}=0.001$ (red
line), $m_{A}/m_{B}=0.01$ (green dashed line) and $m_{A}/m_{B}=0.12$ (blue
dashed-dotted line). The distinct panels correspond to intracomponent
interactions (a) $g_{A}=0$, (b) $g_{A}=0.5$ and (c) $g_{A}=25$. The insets
show the respective pair probabilities. In all cases $g_{AB}=8$. Figure
adapted from (Pflanzer et al., 2009)

###  4.5 Tunneling with mass imbalance

A typical observable for monitoring the tunneling behavior is the population
imbalance among the left (L) and right (R) sites of the double-well in the
course of the evolution. This is easily quantified by the expectation value of
the one-body density given, for example, in the right well (0, $x_{R}$) by

$p_{R}^{\sigma}(t)=\frac{1}{N_{\sigma}}\int_{0}^{x_{R}}dx\rho_{1}^{\sigma}(x,t).$
(32)

This provides a single-particle measure of the tunneling dynamics. The one-
body density of the $\sigma=A,B$ component is given by
$\rho_{1}^{\sigma}(x,t)=\braket{\Psi(t)}{\hat{\Psi}_{\sigma}^{\dagger}(x)\hat{\Psi}_{\sigma}(x)}{\Psi(t)}$,
where $\ket{\Psi(t)}$ is the many-body wave function of the binary system and