$N(\mu)=n_{0}\sqrt{-2/\mu_{0}}[\ln((1+\sqrt{\mu/\mu_{0}})/(\sqrt{1-\mu/\mu_{0}}))-\sqrt{\mu/\mu_{0}}]$.

The solution of Eq. (71) originates from the balance of the cubic self-
repulsion and quadratic attraction of Eq. (70). It occurs within the interval
$\mu_{0}<\mu<0$ and exhibits a flat-top shape for
$\mu-\mu_{0}\ll\left|\mu_{0}\right|$ of size
$L\approx(-2\mu_{0})^{-1/2}\ln[1-(\mu/\mu_{0})^{-1}]$. As a consequence, the
size of the droplet diverges for $\mu=\mu_{0}$ and thus its wave function
tends to a delocalized plane-wave solution of maximally uniform density
$n=n_{0}$. This implies the incompressibility of the BEC which can be
accordingly considered as a fluid. Importantly, the full-width-at-half-maximum
of the quantum droplet
$n(x=0)/2=n_{0}(\mu/\mu_{0})^{2}(1+\sqrt{1-(\mu/\mu_{0})})^{-2}$ decreases for
increasing $\mu$ reaching the value $\sim 2.36/\sqrt{\mu_{0}}$ at $\mu\approx
0.776\mu_{0}$ and afterwards as $\mu\to 0^{-}$ it becomes larger and diverges
as $\sim 1.71/\sqrt{-\mu}$.

Finally, let us note that the above-described exact quantum droplet solution
exists also in the range $\delta g/g<0$ and features a crossover from a
Korteweg-de Vries (Lamb jr, 1980; Hirota and Satsuma, 1981)
$\sim\operatorname{sech}^{2}(\sqrt{-\mu/2x})$ to a bright-like soliton
$\sim\operatorname{sech}(\sqrt{-2\mu}x)$ (Mithun et al., 2020). Other types of
exact solutions of the single-component GGP, corresponding for instance to
dark and anti-dark solitons or kink structures, as well as for the asymmetric
mixture described by Eq. (69) and, for instance, featuring Gaussian-like
shapes or expressed in terms of Jacobi elliptic functions, have been addressed
in some detail in (Mithun et al., 2020). Their stability properties are still
an open issue.

###  8.2 Stationary properties of one-dimensional droplets

The ground state phase diagram in free space and the characteristic properties
of quantum droplets including correlation functions and surface tension have
been studied for a symmetric mixture in (Parisi et al., 2019; Parisi and
Giorgini, 2020) using a numerically exact Diffusion Monte Carlo (DMC) method.
In the large particle number (thermodynamic) limit it was demonstrated that
there is a critical ratio $\left|g_{12}\right|/g\approx 0.45$ above which the
ground state energy of a uniform system exhibits a minimum at a finite
density. This situation corresponds to the liquid phase which can be either
stable or unstable, breaking into droplets above and below the spinodal point,
respectively, where the compressibility diverges. Below this critical value
the energy increases for larger densities, which signals the presence of a gas
phase. Here, the energy minimum corresponds to half of the binding energy of
dimers i.e. $\epsilon_{b}/2=-\hbar^{2}/(2m\tilde{a}^{2})$ with
$\tilde{a}=2\hbar^{2}/(m\left|g_{12}\right|)$ being the 1D scattering length
of the interspecies contact potential and it is associated with a vanishing
density. Note also that this value of $\left|g_{12}\right|/g$ is in accordance
with the prediction of the four-body scattering problem and in particular the
crossing point of the effective interaction between dimers from repulsive to
attractive values (Pricoupenko and Petrov, 2018).

Comparing these exact results with the GGP predictions, an adequate agreement
has been demonstrated in high density regions but large deviations have been
found at small densities where the weak coupling theory is not valid.
Moreover, by inspecting suitable pair correlation functions it was found that
there are interaction regimes of strongly correlated liquid states which lie
beyond the applicability of the GGP framework. Notice also that based on
Bogoliubov theory the phase diagram and the thermodynamic properties of 1D
weakly interacting liquids in terms of the intercomponent interaction and the
temperature have been discussed in (De Rosi et al., 2021). It was shown that
the dominant contribution to the thermal effects stems from the excitation of
the soft density mode.

Figure 36: Ground-state density profiles of 1D droplets for varying particle
number of a symmetric mixture as predicted within the DMC approach (solid
lines). Shown are the cases of $\left|g_{12}\right|/g=0.95$ (left panel) and
$\left|g_{12}\right|/g=0.6$ (right panel). Dashed lines represent the density
obtained in the GGP framework, whilst the dotted line corresponds to the
equilibrium density of a uniform mixture at $\left|g_{12}\right|/g=0.6$ as
calculated using DMC. The normalization parameter ${\rm a}=-2\hbar^{2}/(mg)$
refers to the 1D intracomponent scattering length. Figure from (Parisi and
Giorgini, 2020).

Interestingly, a quantitative comparison between the DMC and the GGP approach
on droplet properties as captured by their density profiles, surface tension
and breathing mode frequency yielded a very good agreement for symmetric
mixtures consisting of a finite particle number (Parisi and Giorgini, 2020).
Deviations between the results of the two methods were mainly identified
within interaction regions where the GGP theory is not expected to be
accurate, i.e. $\left|g_{12}\right|\ll g$, see for instance Fig. 36.

For finite particles numbers, the effect of the surface tension is important
leading to the formation of a self-bound nonuniform configuration. The surface
tension is in general contained in the energy density of the system. In 1D and
for large particle numbers the energy per particle is expressed with respect
to the bulk energy $E_{B}$, which provides the energy per particle of a
homogeneous system in the thermodynamic limit, and the surface energy $E_{S}$
is the first order correction due to the presence of a finite particle number,
leading to

$\frac{E}{N}\approx E_{B}+\frac{E_{S}}{N}+\dots.$ (72)

Higher-order correction terms are neglected here, for simplicity. The surface
tension of the droplet is determined by the ratio of the surface energy
$E_{S}$ to the surface area, which in 1D corresponds to only two points.
Therefore this ratio takes the value $E_{S}/2$.

Figure 37: (a) Droplet ‘wave function’ [Eq. (73)] features a transition from
soliton-like to flat-top configurations when varying $\gamma$. (b) The ratio
of discrete Bogoliubov frequencies $\omega_{\eta}$ to the particle-emission
threshold $-\mu$ with respect to $\gamma$. Red crosses signify the so-called
branching points at which the discrete modes cross $-\mu$ and enter the
continuum of excitations. Dashed and dotted curves represent asymptotic
approximations of the spectrum and the energy of the bound state respectively.
Figure adapted from (Tylutki et al., 2020).

This linear dependence of the energy $E$ on $1/N$ is captured for large $N$
both within the DMC approach and the GGP theory. However, it has been argued
that the GGP underestimates (overestimates) $E$ ($E_{S}$) for
$\left|g_{12}\right|/g<0.95$. Regarding the density of the mixture it has been
demonstrated (Parisi and Giorgini, 2020) that, for a constant ratio
$\left|g_{12}\right|/g$ and increasing particle number, the density of the
mixture becomes higher and wider, while in the large $N$ limit the density
near the center of the cloud becomes flat reaching the equilibrium value of
the uniform system, see Fig. 36. Therefore, an important result is that for
large $N$, the density of the bulk is independent of the particle number
whereas the radius of the cloud continues to expand. The GGP framework shows
differences compared to the DMC predictions already for
$\left|g_{12}\right|/g\approx 0.95$, which are attributed to residual beyond
mean-field effects that are not included in the LHY correction. However, for
$\left|g_{12}\right|/g\ll 1$ these deviations become significant since
interparticle correlations are important and lie far beyond the GGP
description.

Another important characteristic of the droplets is their excitation spectrum,
which has been derived by (Tylutki et al., 2020; Zin et al., 2021; Salasnich,
2018) using a Bogoliubov-de-Gennes linearization analysis. In (Tylutki et al.,
2020) the control parameter is the dimensionless constant

$\gamma=\frac{\pi N^{2/3}\delta g}{2g}.$ (73)

It allows one to regulate the shape of the droplet, see Fig. 37 (a). When
decreasing $\gamma$, the density of the system changes from a flat-top profile
($\gamma\gg 1$) to a Gaussian-like shape as $\gamma$ takes negative values,
gradually approaching a soliton-like configuration for $\gamma\to-\infty$. The
behavior of the mode frequencies as function of $\gamma$ has also been
studied. The conclusion is that the lowest non-trivial collective mode, namely
the breathing mode, remains bounded. It is therefore essentially the only mode
that remains supported by the droplet configuration independently of $\gamma$,
while the higher-lying modes eventually cross into the continuum for
decreasing $\gamma$, see the spectrum in Fig. 37 (b). In the vicinity of the
aforementioned crossings the corresponding mode is related to a large
probability amplitude of finding a particle outside the droplet. More
precisely, flat-top droplets occuring for $\gamma\gg 1$ and having a large
bulk region support plane-wave phonons and resemble constant-density elastic
media with free ends. Decreasing $\gamma$ results in a reduced bulk and the
phonons with the smaller wavelengths approach the continuum faster. At the
balance point, i.e. at $\gamma=0$, the system is described by a GPE with only
a quadratic nonlinearity, which gives rise to a Korteweg-de-Vries type droplet
wave function of the form $\sim\cosh^{-2}(\sqrt{-\mu/2}x)$. A detailed
discussion on the properties of the Korteweg-de-Vries equation and its soliton
solutions can be found in (Hirota and Satsuma, 1981; Hirota, 1972; Wadati and
Toda, 1972; Dodd et al., 1982).

The breathing frequency of the droplet has been estimated using the GGP and
the DMC methods (Parisi and Giorgini, 2020) and significant deviations between
these two approaches have been found only for small particle numbers. Keeping
$\left|g_{12}\right|/g$ fixed, the breathing frequency in the GGP framework
increases for larger particle numbers, reaches a maximum at a certain $N$, and
then decreases approaching zero for very large $N$. As such, it has been
argued that the origin of the breathing mode changes from a surface mode
triggered by surface tension effects to a bulk density mode due to the bulk
compressibility.

Having analyzed the droplet properties in free space a remaining question is
the impact of the external confining geometry. At this point it is worth
mentioning that the effect of an external harmonic trap on the droplet
properties and in general the LHY contribution is still not well understood
(Mistakidis et al., 2021). However, the ground state phase diagram of these
structures for symmetric bosonic mixtures ($N_{A}=N_{B}\equiv N/2$) in 1D
lattice potentials of unit filling $N/2=l$ has been investigated in (Morera et
al., 2020, 2021a) using the DMRG approach in order to simulate the underlying
Bose-Hubbard Hamiltonian, see also Section 9 and in particular the discussion
around Eq. (93). The tight-binding Bose-Hubbard model that describes the
system under consideration is given by

$\begin{split}H=-t\sum_{i}\sum_{\alpha=A,B}(\hat{b}_{i,\alpha}^{\dagger}\hat{b}_{i+1,\alpha}+\hat{b}_{i,\alpha}\hat{b}_{i+1,\alpha}^{\dagger})+\frac{U}{2}\sum_{i}\sum_{\alpha=A,B}[\hat{n}_{i,\alpha}(\hat{n}_{i,\alpha}-1)]+U_{AB}\sum_{i}\hat{n}_{iA}\hat{n}_{iB},\end{split}$
(74)

where $\hat{n}_{i\alpha}=\hat{b}_{i\alpha}^{\dagger}\hat{b}_{i\alpha}$ is the
number operator of species $\alpha=A,B$ acting at lattice site $i$. The
corresponding bosonic creation and annihilation operators are denoted by
$\hat{b}_{i\alpha}^{\dagger}$ and $\hat{b}_{i\alpha}$. For simplicity reasons,
the same tunneling strength $t>0$ and repulsive intraspecies interaction $U>0$
is considered for both species, with $t$ being kept fixed to define an energy
scale for the system. The attractive interspecies interaction is described by
$U_{AB}<0$, obeying $\left|U_{AB}\right|<U$. Moreover, for the
characterization of the emergent droplet properties the dimensionless ratio
$r=1+U_{AB}/U>0$ is introduced such that the balance point of the mean-field
interactions occurs at $r=0$.

Figure 38: (Left panel) Phase diagram of lattice trapped 1D droplets as
predicted by the two-component Bose-Hubbard model of Eq. (74) for
$N_{A}=N_{B}=L=64$ using DMRG. Green triangles mark the MI-PSF phase
transition, the dashed curve denotes the MI-PSF transition line, the orange
triangles show the PSF-2SF transition, while blue squares indicate the
parameter region where droplets having a smaller size than the lattice length
occur. (Right panels) Density profiles in the weakly interacting regime where
$r=0.01$ for different particle numbers. A gaseous 2SF state characterized by
$U/t=0.1$ (right top panel) is presented and a D-2SF (right bottom panel) with
$U/t=2.5$ respectively. The black dashed lines represent fits using Eq. (75).
In the lowest panel an asymmetric configuration with $U_{AA}=5U_{BB}$ is
depicted using the red crosses. Figure adapted from (Morera et al., 2020).

The phase diagram of the 1D droplet setting for varying interactions is
presented in the top panel of Fig. 38. For strong couplings $U\gg t$ the Mott-
insulator (MI) and the pair superfluid (PSF) phases emerge in the system.
Their transition point takes place at $(r)_{c}\approx 9.25t^{2}/U^{2}$ and it
is identified by the behavior of the related pair correlation function which
decays exponentially in the MI and algebraically in the PSF. Note that in the
MI state coherence is exponentially lost (Greiner et al., 2002; Stöferle et
al., 2004). However, the PSF is related to the formation of atom pairs
exhibiting long-range phase coherence and in particular a power-law decay of
the off-diagonal quasi-long-range order222222Recall that 1D quantum phases
cannot exhibit true long-range order, see also Section 6.4. Rather quasi-long-
range order appears meaning that the correlation function
$\braket{\Lambda^{\dagger}(x)\Lambda(0)}$ of an operator $\Lambda(x)$ falls of
algebraically $\sim\left|x\right|^{\beta}$ as $\left|x\right|\to\infty$. (Hu
et al., 2009; Kuklov et al., 2004b, a), absence of phase coherence in the one-
body correlator and the appearance of a finite energy gap associated with spin
excitations.

In order to capture the self-bound configurations, one can consider the
density profile of the system and use a fit to the following symmetrized Fermi
function 232323The symmetrized Fermi function of Eq. (75) is constructed in
order to avoid the cusp of the well-known Fermi function at $x=0$. Relevant
elaborated discussions regarding its properties can be found e.g. in Refs.
Sprung and Martorell (1997); Buck and Pilt (1977); Grammaticos (1982).

$n_{i}=\frac{n_{M}\sinh[R/(2a)]}{\cosh[R/(2a)]+\cosh(i/a)},$ (75)

where $n_{M}$, $R$ and $a$ are free parameters. Evidently, this function has a
flat-top profile of size $R$ and shows an exponential decay at its edges of
typical length scale $a$. This is shown in the right hand side panel in Fig.
38. The size of the observed droplets for fixed $r$ depends strongly on $U$.
Particularly, for a specific $r$ there exists a value of $U/t$ above which the
system resides in a MI state extending over the entire lattice. However,
decreasing $U/t$, the PSF phase is entered where the droplet size shrinks up
to a critical value of $U/t$. For weaker interactions a two-atomic superfluid
(2SF) state is realized, which is characterized by an increasing droplet size
that becomes comparable to the system size in the noninteracting limit. Within
the 2SF state each species exhibits separately quasi-long phase coherence.
Therefore, the transition MI-PSF is identified by the appearance of quasi-
long-range coherence in the two-body correlators while the PSF-2SF phase
transition is captured by the build-up of quasi-long-range coherence in the
one- and two-body correlators. Note that the PSF-2SF transition can occur
either when both phases are in their gaseous (D-PSF) or droplet (D-2SF)
configuration, where in the latter case the density profile has a smaller
extension compared to the former. The impact of the filling factor of the
lattice on the density profiles of the phases is discussed in (Morera et al.,
2020), but a detailed analysis of this issue is still an outstanding question.
Calculations of the bipartite entanglement entropy of the droplet
configurations have revealed that for increasing particle number the
entanglement becomes maximal when the droplet size [see also Eq. (75) and left
panels of Fig. 38] is at a minimum, while it decreases for a larger number of
particles. As such, droplets with a minimum size exhibit more pronounced
quantum effects.

The structure of the lattice trapped mixture described by Eq. (74) in the
vicinity of the mean-field balance point $r=1+U_{AB}/U\ll 1$ with strong
interactions $\left|U_{AB}/t\right|\approx U/t\gg 1$ was addressed in (Morera
et al., 2021a). For the description of this system an effective dimer model
was used (Kuklov et al., 2004a, b):

$H_{\text{eff}}^{D}=-J^{(2)}\sum_{n}(\hat{c_{n}}^{\dagger}\hat{c}_{n+1}+\hat{c_{n}}\hat{c}_{n+1}^{\dagger})+\frac{U^{(2)}}{2}\sum_{n}\hat{N}_{n}^{D}\hat{N}_{n}^{D}+V^{(2)}\sum_{n}\hat{N}_{n}^{D}\hat{N}_{n+1}^{D},$
(76)

where $\hat{N}_{n}^{D}$ is the dimer number operator and $\hat{c}_{n}$
($c_{n}^{\dagger}$) is the dimer annihilation (creation) operator obeying
$[\hat{c}_{n},\hat{c}_{m}^{\dagger}]=\delta_{n,m}$. The parameter
$J^{(2)}=2t^{2}(1+r)/U$ describes the hopping of dimers, while
$U^{(2)}=Ur-4t^{2}r/U$ and $V^{(2)}=-4t^{2}(1-r)/U$ correspond to the
interaction strength between two dimers and the nearest-neighbor interaction,
respectively.

From $H_{\text{eff}}^{D}$ analytical estimates of the dimer-dimer effective
range and scattering length, as well as the tetramer bound state threshold and
binding energy can be derived (Pricoupenko and Petrov, 2018; Hadizadeh et al.,
2011; Greene et al., 2017). The $s$-wave scattering length $a_{DD}$ of the
two-body problem of dimers can be expressed as (Valiente and Petrosyan, 2009)

$\frac{a_{DD}}{d}=\frac{(r-1)rU^{2}/t^{2}-4(r^{3}+3r+4)}{2r^{2}U^{2}/t^{2}-8},$
(77)

where $d$ is the lattice constant. From the formula one can see that when
$r_{c}=2t/U$ the scattering length diverges. Note that it also changes its
sign from $a_{DD}\to+\infty$ to $a_{DD}\to-\infty$ with increasing $U/t$ while
keeping $r$ constant. The effective dimer-dimer coupling constant of the form
$U_{DD}=-2\hbar^{2}/(m_{D}^{*}a_{DD})$ tends to zero at $r_{c}$ while it is
attractive (repulsive) for $a_{DD}>0$ ($a_{DD}<0$). Here $m_{D}^{*}$ denotes
the effective mass of the dimer. Moreover, in the attractive regime a dimer-
dimer bound state, the so-called tetramer, forms in the vicinity of $U_{DD}=0$
with a binding energy given by

$E_{B}\approx-\frac{\hbar^{2}}{m_{D}^{*}a_{DD}^{2}}.$ (78)

The predictions of the effective model from Eq. (76) were compared with DMRG
calculations of the full Hamiltonian and were found to be in a very good
agreement especially in the deep dimer regime characterized by $U/t\gg 1$, see
also Ref. (Valiente and Petrosyan, 2009; Pricoupenko and Petrov, 2018).
However, for $U/t\ll 1$, the effective model overestimates the depth of the
bound state in contrast to the DMRG results. Remarkably, the transition from a
gas ($U_{DD}>0$) to a liquid ($U_{DD}<0$) phase has been found in the
interaction region where the formation of a bound tetramer occurs according to
Eq. (78), i.e. when the effective dimer-dimer interactions change from
repulsive (gas) to attractive (liquid). Thus, the transition threshold
exhibits a universal behavior with respect to the relevant parameters of the
dimer-dimer scattering problem (Pricoupenko and Petrov, 2018), i.e. the
scattering length and the effective range. The stability of the liquid phase
is related to the nonzero effective range of the dimer-dimer scattering
problem $r_{e}/d\approx(1-r_{c})(3r_{c}+1)/(r_{c}+1)^{2}$. This result is in
sharp contrast to the un-trapped scenario where the underlying stabilization
mechanism is caused by the three-dimer interactions (Guijarro et al., 2018).
Finally, entering the regime of sufficiently strong dimer-dimer repulsions
($a_{DD}\ll d$) the fermionization of the dimers forming an effective Tonks-
Girardeau state can be seen.

###  8.3 Dynamics of one-dimensional droplets

To further understand droplet properties and their potential applications, one
studies dynamics, e.g. by considering collisions between two droplets
(Astrakharchik and Malomed, 2018) or their spontaneous generation trigerred by
an instability mechanism (Mithun et al., 2020). We should note here that the
relevant dynamical considerations are restricted to only a few droplets and
thus the current understanding is still far from complete.

The collisional dynamics between two counterpropagating droplets is very
interesting since it allows one to distinguish them from bright soliton
structures in spite of the fact that droplets have density profiles similar to
those of bright solitons. Recall that a soliton maintains its shape when
moving with constant velocity and also when colliding with another soliton as
long as the system remains nearly integrable (Kevrekidis et al., 2007).
Similarly, droplets whose wave function is described by Eq. (71), exhibit a
constant density profile $n(x;t)=\left|\Psi(x;t)\right|^{2}$ in the course of
their time-evolution if moving with a constant velocity. However, they
generally do not preserve their shape when undergoing pairwise collisions
(Astrakharchik and Malomed, 2018), a result that clearly differentiates them
from bright solitons. Let us also remark here another interesting property
that concerns the stability of a single droplet when it is exposed to a
periodic density modulation $\sim\cos(kx)$ of wavenumber $k$. It has been
shown that there is a particle dependent critical wavenumber $k_{c}$ above
which the droplet splits into several pieces propagating in opposite
directions while for $k<k_{c}$ it maintains its shape featuring almost
periodic oscillations.

Figure 39: Density plots of two colliding droplets whose initial form is given
by Eq. (79). Different situations for symmetric ($N_{1}=N_{2}$) and asymmetric
($N_{1}\neq N_{2}$) mixtures are shown for different incident momenta $k$ (see
legends). Unless stated otherwise, the collisions are in-phase ($\phi=0$).
Figure from (Astrakharchik and Malomed, 2018).

Two counterpropagating droplets were considered in (Astrakharchik and Malomed,
2018) assuming that an initial ‘wave function’ has the form

$\Psi(x;0)=e^{ikx/2+\phi}\Psi_{1}(x+x_{0})+e^{-ikx/2}\Psi_{2}(x-x_{0}),$ (79)

where $\Psi_{1}(x)$, $\Psi_{2}(x)$ are given by Eq. (71), $\pm x_{0}$ denote
the initial locations of the droplets, $\pm k/2$ corresponds to their initial
momenta and $\phi$ is the relative phase between them. (Astrakharchik and
Malomed, 2018) demonstrated that the collisional properties of the droplets
depend strongly on the considered particle number, their velocity and relative
phase. As discussed in Section 8.2, for small particle numbers the droplet
profile has a Gaussian-like shape while for large atom numbers it is of a
“puddle” type. Below, the basic collisional features reported in
(Astrakharchik and Malomed, 2018) are summarized, see also Fig. 39.

Two small droplets preserve their shapes after the collision and show an
interference pattern at the collision point, which is similar to the interplay
of coherent matter waves. In sharp contrast, two counterpropagating “puddles”
generally display aspects of inelastic collisions, accompanied by mass-
transfer phenomena and leading to a larger number of droplets or a merger of
the outgoing ones for fast and slow velocities, respectively. Depending on the
phase difference it is possible to tune the degree of excitation and the
population of the outcoming droplets, especially for initially large droplets.
On the other hand, the collision of a “puddle” and a Gaussian-like droplet
results in the excitation of the former performing internal vibrations and an
almost un-perturbed shape of the latter. Two particle imbalanced Gaussian-like
droplets show an almost elastic collision with their outgoing trajectories
being only slightly affected.

The impact of beyond LHY correlations in the droplet nonequilibrium dynamics
induced by an interspecies interaction quench within the nonperturbative ML-
MCTDHX approach was studied in (Mistakidis et al., 2021). Specifically, a
sudden change from weak to strong intercomponent attractions enforces the
droplet expansion acquiring flat top signatures during the dynamics. Here, the
expansion velocity is found to be reduced as compared to the GGP prediction.
Following the opposite quench, namely from strong to weak interactions, a
breathing motion is triggered whose frequency becomes slightly smaller than in
the GGP framework due to residual higher-order correlations (Parisi and
Giorgini, 2020; Astrakharchik and Malomed, 2018). In this sense, both the
expansion and the breathing dynamics of the droplet provide experimental means
for identifying beyond LHY physics. Interestingly, two-body correlations
build-up in the course of the evolution evidencing that two bosons of the same
component experience a bunching at long distances and an anti-bunching at the
same location.

Figure 40: Time-evolution of $\left|\Psi(x,t)\right|$ representing an
initially perturbed plane-wave solution in the ultracold regime (left panel)
and for the case of the presence of dissipation (condensate thermal fraction)
with strength $\gamma=0.03$ (right panel). The spontaneous generation of an
array of quantum droplets and their subsequent coalescence can be observed.
When dissipation is taken into account to emulate thermal effects of the
background its fluctuations are minimized. The initial homogeneous density
corresponds to $n\approx 0.1$ while $\delta g/g=0.05$. Figure adapted from
(Mithun et al., 2021).

The nucleation of quantum droplets from plane-wave solutions due to the
manifestation of the modulation instability242424The modulational instability
wave phenomenon is related with the exponential growth of perturbations due to
the presence of nonlinearity and anomalous dispersion. It manifests in various
physics areas e.g. ranging from BECs (quantum fluids) (Nguyen et al., 2017;
Everitt et al., 2017), and waveguides (Tai et al., 1986) to systems in
nonlinear optics (Agrawal, 2000) and charged plasmas (Thornhill and Ter Haar,
1978). has been discussed in (Mithun et al., 2020, 2021). In particular, let
us consider a stationary plane-wave subjected to small perturbations i.e.
$\Psi(x;t)=[\sqrt{n}+\delta\Psi(x;t)]e^{-i\mu t}$. Here,
$\delta\Psi=\zeta\cos(kx-\Omega t)+i\eta\sin(kx-\Omega t)$ with $\zeta$,
$\eta$ being infinitesimal amplitudes, while $k$ and $\Omega$ denote the
wavenumber and frequency of the pertubation. Performing a linearization
analysis with respect to $\delta\Psi$ yields the corresponding Bogoliubov-de-
Gennes equation of motion and leads to the dispersion relation

$\Omega^{2}=\frac{k^{4}}{4}+\bigg{(}\frac{\delta
g}{g}n-\frac{\sqrt{n}}{\sqrt{2}\pi}\bigg{)}k^{2}.$ (80)

The modulation instability occurs when $\Omega$ becomes imaginary, namely when
$n<[2\pi^{2}(\delta/g)^{2}]^{-1}$ or $k^{2}<4[\sqrt{n}/(\sqrt{2}\pi)-(\delta
g/g)n]\equiv k_{0}^{2}$. This imaginary contribution becomes maximal for
$k_{max}=k_{0}/\sqrt{2}$. Letting the perturbed plane-wave solution evolve
within the GGP framework of Eq. (70) and operating within parameter regions
where the modulation instability is expected to take place [see Eq. (80)]
leads to the emergence of chains of quantum droplets, see the left hand side
panel of Fig. 40. Their number reduces during the evolution due to mergers of
the colliding droplets and these observations clearly indicate that the
modulation instability is indeed the underlying mechanism for the creation of
quantum droplets in this setting. This result also holds true for asymmetric
mixtures characterized either by a population imbalance or by different
intracomponent coupling strengths. Finally, in an attempt to account for
thermal effects, a phenomenological dissipative term of strength $\gamma$ was
considered in the GGP by (Mithun et al., 2021). It was shown that in the
presence of dissipation the overall generation of droplets is similar and
happens on the same timescale as the zero temperature scenario. However, in
this situation the formation of "large" droplets becomes more prominent and
the background fluctuations are minimized, see the right hand side panel of
Fig. 40. It should be emphasized at this point that the role of thermal
fluctuations and their adequate modelling including their impact on the LHY
term is certainly an issue that deserves further study.

Independently, the quench dynamics of a quasi-1D repulsively interacting
binary bosonic mixture was investigated in (Cidrim et al., 2021). The quench
protocol consisted of a linear ramp of the interactions over a finite time
interval and the calculations were based on the GGP framework including also
the effect of three-body losses (Semeghini et al., 2018). It was found that
after quenching to weak (strong) intercomponent couplings the arising
modulation instability allows for the formation of dark soliton trains (self-
bound droplets). An analytical prediction of the number of generated solitons
was also provided through the unstable eigenmode of the Bogoliubov spectrum.
In sharp contrast, a quench from the weakly attractive (soliton) regime to
strong interactions (droplet region), produces a splitting of the original gas
into two bright solitons, but not a train (see also (Cheiney et al., 2018)).

Over the past few years it has been demonstrated that quantum droplets feature
inelastic collisions and that the modulation instability process provides
means for their spontaneous generation. Despite the fact that various
theoretical questions for the dynamical response of droplet configurations in
1D still remain open, it is worthwhile to comment that these systems
constitute fruitful platforms for probing beyond mean-field effects in weakly
interacting ultradilute Bose gases (Navon et al., 2011). They also provide a
new arena for generalizing already well studied concepts such as non-linear
excitations and impurity physics. Another future research direction of
particular importance is the study of the competition between temperature
effects and quantum fluctuations and also the generic modelling of temperature
(Mithun et al., 2021; De Rosi et al., 2021). Finally, we note that most of the
current investigations have focused on the dynamics of quantum droplets
described by the GGP framework. However, beyond LHY correlations result in a
reduced velocity of expanding droplets and a slightly smaller breathing
frequency. This motivates beyond LHY studies to reveal higher-order
correlation effects in the course of the time evolution.

##  9 Methods: review

###  9.1 Exact Diagonalization

Arguably, the most widespread numerical technique to approach quantum systems
with a few degrees of freedom is exact diagonalization (ED). ED is a term that
applies broadly to any numerical method based on writing the Hamiltonian
matrix in a prescribed basis and then diagonalizing the obtained matrix. The
ED method is useful if the Hilbert space is sufficiently252525The meaning of
the word ‘sufficiently’ here strongly depends on the problem at hand, in
particular, on the sparsity of the Hamiltonian matrix. small. It is widely
used to analyze lattice many-body problems (see, e.g., (Weiße and Fehske,
2008)), and few-body systems. The main advantage of ED is the possibility to
control the accuracy of calculations – one can reduce the error bars in a
systematic manner. In addition, one can obtain the full energy spectrum of the
Hamiltonian, which gives a basis for simulating time dynamics, see, e.g.,
recent studies on time dynamics of few-body systems (Ebert et al., 2016; Yin
et al., 2016; Nandy and Sowiński, 2020, 2021).

For lattice systems with a finite number of degrees of freedom the Hamiltonian
matrix is in principle finite. However, the ‘exponential curse’, i.e., the
fact that a Hilbert space will generally grow exponentially with the number of
particles, means that in practice it is always necessary to perform some kind
of truncation procedure to keep things manageable. A number of general purpose
codes exist in the open source domain, including the widely used ALPS (Bauer
et al., 2011) and QuSpin packages (Weinberg and Bukov, 2017, 2019).

For continuum systems, which are the focus of this review, the Hilbert space
is infinite to begin with, and, to employ ED, one works with a truncated
basis. One common strategy is to use a non-interacting one-body basis with
$n_{\mathrm{max}}$ states. To ensure a converged result when working with
finite bases, it is customary to calculate the spectra, eigenstates, and other
properties of the system as a function of the basis size. A key question is
the approach to the infinite basis size limit, i.e., how fast will the results
converge as the basis size is changed.

For a concrete setting with a finite number of particles there are a number of
systematic tricks that can be applied to either speed up convergence or obtain
higher accuracy for a given quantity of interest, or both. This is a vast
field that features a lot of creativity in ‘clever’ basis choices that
accomplish these tasks more efficiently. We will not attempt a full review,
but merely point out some results that have been successful in cold atoms, and
in particular in the direction of few-body static and dynamic studies. The
focus here is on short-range interactions.

For particles interacting via zero-range potentials in one spatial dimension,
it is known that convergence as a function of $n_{\mathrm{max}}$ can be slow.
The slow convergence is due to a cusp of the wave function. It leads to a non-
analytical behavior, which is hard to reproduce using bases based upon non-
interacting models. For example, for eigenfunctions of harmonic oscillator
potential, the convergence is expected to be $1/\sqrt{n_{\mathrm{max}}}$
((Grining et al., 2015)), whereas for a plane wave basis the convergence
should be $1/n_{\mathrm{max}}$ ((Volosniev and Hammer, 2017b; Jeszenszki et
al., 2018)). The difference between $1/\sqrt{n_{\mathrm{max}}}$ and
$1/n_{\mathrm{max}}$ is related to a different dependence of the maximal one-
body energies on $n_{\mathrm{max}}$, which in those cases determine the
resolution of non-interacting wave functions. It is worth noting that for
attractive interactions, the slow convergence pushes the solution out of reach
faster than for repulsive interactions (see, e.g., (D’Amico and Rontani,
2014)).

To improve convergence, one can introduce an effective interaction instead of
the bare one. Ideally, the effective interaction is an interaction term that
produces results close to the exact ones already for relatively small values
of $n_{\mathrm{max}}$. One simple empirical method to obtain an effective
interaction from the bare one is to rescale the interaction strength as

$g\to\frac{g}{1+g/g_{0}},$ (81)

where $g_{0}$ is an empirical constant that depends on $n_{\mathrm{max}}$,
such that $g_{0}\to\infty$ as $n_{\mathrm{max}}\to\infty$. It was shown that a
rescaling of Eq. (81) significantly improves the accuracy of the energies
(Ernst et al., 2011). Later, (Jeszenszki et al., 2018, 2019) studied
convergence patterns of the approach, and showed that renormalization allows
one to improve convergence by an order of magnitude, and also provided a
transformation method to improve the wave functions around cusps for strong
interactions (Jeszenszki et al., 2020).

Another effective interaction for cold atoms can be obtained using the Lee-
Suzuki method from nuclear physics (Christensson et al., 2009; Rotureau,
2013). For studies of few-body systems in one spatial dimensions, see
(Lindgren et al., 2014; Dehkharghani et al., 2015a; Volosniev and Hammer,
2017b); freely available software is introduced in (Rammelmüller et al.,
2022). The effective interaction in those works is based upon the solution of
a two-body problem. In a truncated Hilbert space, the effective interaction
can be written as

$V_{\mathrm{eff}}=Q^{T}EQ-H_{0},$ (82)

where $E$ is the matrix that contains two-body energies on a diagonal, $Q$ is
an orthogonal matrix, and $H_{0}$ is a non-interacting two-body Hamiltonian.
The $Q$ is chosen as $u/\sqrt{u^{T}u}$, where $u$ is a matrix whose rows
contain a projection of eigenvectors on the truncated Hilbert space. The two-
body solutions for a ring and harmonic potentials are discussed in Secs. 2 and
3. Using those solutions one can construct an effective potential that
incorporates correct two-body information in a truncated Hilbert space, thus
speeding up the convergence.

Finally, we mention that the convergence can be improved even for a
Hamiltonian with a bare interaction if one modifies the basis used for
diagonalization. For example, one can use a basis with additional variational
parameters, e.g., the frequency of the harmonic potential (Kościk, 2018); a
basis that utilizes conserved quantities (Harshman, 2012); a basis based upon
random states (explicitly correlated Gaussians) that are not related to a non-
interacting problem (Andersen et al., 2016), see also (Mitroy et al., 2013)
for a review of that approach. Those methods can be very successful. However,
they might require a system-dependent pre-selection of suitable basis states
and departure from the established numerical techniques, which stall their
development.

###  9.2 Bethe Ansatz-Based Methods

In Section 2 we presented the form of the Bethe ansatz wave function for a
bosonic system and discussed the solutions in the few-body limit. For general
$N$, assuming a system described by Eq. (3) and imposing the contact
interaction condition on the derivatives, as well as periodic boundary
conditions, yields a set of Bethe ansatz equations (BAE) given by

$k_{j}L=2\pi
I_{j}-(N-1)\pi-2\sum_{l=1}^{N}\text{arctan}\left(\frac{k_{j}-k_{l}}{c}\right),$
(83)

where $I_{i}$ is an integer or half-integer for a total odd or even number of
particles, respectively. Solving these equations specifies all the possible
values of $k_{i}$, thus determining the complete wave function. To express a
physical solution, these values must all be different from each other. It is
worth pointing out here that the solutions strongly depend on the sign of the
interaction: for the repulsive case, all solutions are real (Yang and Yang,
1969), while for the attractive case complex solutions called strings are
allowed (Takahashi, 2005). Solving the eigenvalue equation $H\psi=E\psi$ then
returns the total energy of the system, which is simply given by

$E=\sum_{i}k_{i}^{2}.$ (84)

The problem of finding the whole spectrum of the model is thus reduced to
solving the BAE, from which the ground state energy, elementary excitations
and correlation functions can be derived. Remarkably, the method allows one to
handle both the repulsive and attractive regimes as well as to deal with an
arbitrary number of particles, and has been successfully employed in a number
of applications. A comparison between experimental results obtained for the
pair correlations of a one-dimensional Bose gas and the results predicted by
the Bethe ansatz approach is shown in Fig. 41, where a decrease of pair
correlations can be observed as the interactions are made stronger.

Figure 41: Experimental measurement of the pair correlations in the repulsive
one-dimensional Bose gas. The solid curve shows the theoretical prediction
obtained using the Bethe ansatz for the Lieb-Liniger model (Gangardt and
Shlyapnikov, 2003). Figure from (Kinoshita et al., 2005).

For the fermionic case, the situation becomes more involved since one has to
deal with different component or spins. Thus, the spatial wave function,
provided by the coordiante Bethe ansatz, must be combined with a corresponding
spin wave function. Yang’s solution for the spin wave function in a two-
component system of $N$ fermions with $N_{\downarrow}$ down spins is obtained
with the generalized or Bethe-Yang ansatz, and reads (Yang, 1967; Takahashi,
2005)

$\Phi=\sum_{R}A(R)F_{P}(\Lambda_{R1},y_{1})...F_{P}(\Lambda_{RN_{\downarrow}},y_{N_{\downarrow}}),$
(85)

where the $y_{i}$ are the coordinates of spin-down particles with
$i\in\\{1,...,N_{\downarrow}\\}$ and $R$ are permutations acting on these
coordinates. For each permutation, the functions $F(\Lambda,y)$ depend on the
spin rapidities $\Lambda$ and are given by

$F(\Lambda,y)=\prod_{j=1}^{y-1}\frac{k_{j}-\Lambda+ic^{\prime}}{k_{j+1}-\Lambda-
ic^{\prime}}$ (86)

where $c^{\prime}=c/2$. In analogy with the bosonic case, the Yang-Baxter
equation allows for the solution of the fermionic problem. The periodic
boundary conditions are now satisfied with a set of $N$ equations for the
quasimomenta $k$ and rapidities $\Lambda$

$\displaystyle k_{j}L=2\pi
I_{j}-2\sum_{\alpha=1}^{N_{\downarrow}}\arctan\left(\frac{k_{j}-\Lambda_{\alpha}}{c^{\prime}}\right),$
(87)

and

$\sum_{j=1}^{N}2\arctan\left(\frac{\Lambda_{\alpha}-k_{j}}{c^{\prime}}\right)=2\pi
J_{\alpha}+2\sum_{\beta=1}^{N_{\downarrow}}\arctan\left(\frac{\Lambda_{\alpha}-\Lambda_{\beta}}{c}\right).$
(88)

In these expressions, $I_{j}$ and $J_{\alpha}$ are integer or half-odd
integers, depending on the numbers $N_{\uparrow}$ and $N_{\downarrow}$ in each
component. The energy expression is the same as Eq. (84). However, we now have
a nested set of the Bethe ansatz equations due to the spin degree of freedom.
The spin wave function obtained with the Bethe-Yang ansatz yields the solution
of the two-component fermionic problem when coupled with the spatial wave
function. It is important to point out that these parts must be combined by
taking into consideration the appropriate permutation symmetries, which can be
obtained from the Young tableaux (Guan et al., 2013). The final wave function
must then by antisymmetric upon simultaneous permutation of the coordinates
and spins of two particles.

The Bethe ansatz approach thus offers an elegant way to solve a number of one-
dimensional models of particles with short range interactions. In addition, it
also allows for the solution and discussion of physical properties of other
types of models in the context of ultra-cold atoms - see, for example, (Links
et al., 2003; Santos et al., 2006; Rubeni et al., 2012, 2017). However, the
wave functions produced can often become too complex and of little practical
use in the calculation of experimentally relevant quantities. This is due to
the fact that they grow factorially with the number of particles, making them
hard to handle for increasingly large $N$. This characteristic is particularly
problematic in a dynamical setting, where the effect of integrability is
expected to be more pronounced.

To address these issues, efficient numerical methods based on – or employing
technical aspects related to – the Bethe Ansatz have been developed and
applied to the solutions of integrable many-body problems in recent years. The
first of these, known as ABACUS (Algebraic Bethe Ansatz-based Computation of
Universal Structure factors) (Caux, 2009) is based on utilizing the solutions
provided by the Bethe ansatz, along with an optimized scan of the Hilbert
space to specifically tackle the problem of calculating the generalized
dynamical structure factor, defined as

$S^{a\bar{a}}(k,\omega)=\frac{1}{N}\sum_{j,j^{\prime}=1}^{N}e^{-ik(j-j^{\prime})}\int_{-\infty}^{+\infty}dte^{i\omega
t}\langle\mathcal{O}^{a}_{j}(t)\mathcal{O}^{\bar{a}}_{j^{\prime}}(0)\rangle$
(89)

where
$\langle\mathcal{O}^{a}_{j}(t)\mathcal{O}^{\bar{a}}_{j^{\prime}}(0)\rangle$ is
the ground state expectation value of a given correlation function at sites
$j,j^{\prime}$ and $a,\bar{a}$ denote site-specific indices. Using the Lehmann
spectral representation (Fetter and Rokhsar, 1998), the structure factor can
be cast as

$S^{a\bar{a}}(k,\omega)\frac{2\pi}{N}\sum_{\mu}|\langle\lambda^{0}|\mathcal{O}^{a}_{k}|\mu\rangle|^{2}\delta(\omega-
E_{\mu}-E_{0}),$ (90)

where $|\lambda^{0}\rangle$ and $E_{0}$ denote the ground state and its
correspondent energy. The essential step in obtaining the dynamical structure
factor thus becomes calculating the form factors
$\langle\lambda^{0}|\mathcal{O}^{a}_{k}|\mu\rangle$ and summing over the
intermediate eigenstates $|\mu\rangle$. The form factors can be expressed
analytically using the algebraic Bethe ansatz for different models, such as
the one-dimensional Bose gas (Slavnov, 1989, 1990) and the Heisenberg chain
(Kitanine et al., 1999, 2000). The summation over intermediate eigenstates is
performed by assuming certain properties of the form factors for large systems
and realizing an optimized scan of the Hilbert space. Once the method returns
a list of form factors, a dynamical correlation function can be calculated.

Since its development, ABACUS has been used to obtain accurate results for the
correlation functions of different one-dimensional integrable many-body
models. Applications include the calculation of properties of the one-
dimensional Bose gas such as one-body dynamical correlations, (Caux et al.,
2007), correlations of highly excited (Fokkema et al., 2014), and finite-
temperature regimes (Panfil and Caux, 2014), as well as the four-spinon
structure factor of the Heisenberg model (Caux and Hagemans, 2006).
Additionally, the method has been successfully applied in the prediction of
experimental results such as the dynamical structure factor (Fabbri et al.,
2015) and the excitation spectrum of the one-dimensional Bose gas (Meinert et
al., 2015) (see Fig. 42).

Figure 42: Experimentally measured Bragg-excitation spectra of a one-
dimensional Bose gas at finite temperatures. Panels (a)-(d) display the
results for increasingly stronger interactions. The solid curves show the
averaged theoretical results obtained from the dynamical structure factor
calculated with ABACUS. Figure adapted from (Meinert et al., 2015)

.

A notable limitation of the Bethe ansatz is the fact that it is restricted to
homogeneous systems, where the integrability condition is preserved. Such a
feature can be troublesome because most experimental results with cold atoms
depend on the presence of a (usually harmonic) trapping potential. While the
local density approximation (LDA) is able to reproduce equilibrium results in
the presence of external potentials, that is not the case for dynamics; one
particular result that posed a challenge to the theoretical developments of
its time is the experiment known as Quantum Newton’s Cradle (Kinoshita et al.,
2006). In this experiment, a many-body Bose gas is initialized in an out-of-
equilibrium state, which evolves in time with the particles oscillating and
colliding in a trap. Interestingly, the system does not seem to thermalize
even after several collisions, meaning that the particles are only allowed to
exchange momentum and their momentum distributions do not change in time
considerably. Such a feature is a characteristic of integrable systems, where
diffraction is absent. However, while usual Bethe ansatz-based techniques are
not able to incorporate the presence of the trapping potential in a real-time
simulation of dynamics, the standard hydrodynamic approaches could not
correctly reproduce the lack of thermalization observed in the experiments.

These limitations were overcome with the development of a new tool that became
known as Generalized Hydrodynamics (GHD). It is often described as an emergent
hydrodynamic theory used to obtain the large scale dynamics of integrable
systems far from equilibrium, and was originally applied to the study of
transport properties of the XXZ chain (Bertini et al., 2016), as well as the
sinh-Gordon and Lieb-Liniger models (Castro-Alvaredo et al., 2016). The method
is built upon principles such as the generalized Gibbs ensemble (Rigol et al.,
2007), which describes the dynamics and relaxation of integrable systems to
steady states after a long time evolution, and local entropy maximization,
defined utilizing the language of the thermodynamic Bethe ansatz. The main
goal of GHD is to solve the Euler-scale equations

$\partial_{t}\rho_{p}(p,x,t)+\partial_{x}\left[v^{\text{eff}}(p,x,t)\rho_{p}(p,x,t)\right]=0$
(91)

where $\rho_{p}(p,x,t)$ denotes neither a fluid density in the classical sense
nor a density of particles; instead, it describes a density of
“quasiparticles” moving with an effective velocity $v^{\text{eff}}$. This
density thus defines the state of a fluid cell, which depends on momentum,
time and position. Since its original proposal, GHD has been developed in
great detail (De Nardis et al., 2018; Bertini et al., 2018; Alba et al., 2019;
Nardis et al., 2019; Vu and Yoshimura, 2019; Bastianello et al., 2019; Bertini
et al., 2021) and applied to a number of different many-body models (Doyon et
al., 2018; Bastianello et al., 2018; Doyon, 2019; Urichuk et al., 2019;
Mestyán et al., 2019; Nozawa and Tsunetsugu, 2020; Bettelheim, 2020),
including the Quantum Newton’s Cradle problem (Caux et al., 2019). The method
has been experimentally verified with measurements of the time-evolution of
the density profile of an Rubidium gas trapped in an atom chip (Schemmer et
al., 2019) (see Fig. 43).

Figure 43: Experimental measurement of the expansion of an atomic cloud of
87Rb released from a double-well potential generated on an atom chip. The
solid curves show the GHD results. Figure from (Schemmer et al., 2019).

###  9.3 Density Matrix Renormalization Group

Very few numerical methods for quantum systems are able to address static and
dynamic problems outside the realm of integrable models while still being
capable of handling fairly large systems. A remarkable exception is the case
of the Density Matrix Renormalization Group (DMRG) method (Schollwöck, 2005,
2011), originally proposed by White in (White, 1992, 1993). It quickly proved
to be one of the most powerful and versatile methods for obtaining the low-
energy properties of one-dimensional quantum many-body models defined on
lattices.

DMRG is an iterative and variational numerical method. Its core idea is
relatively simple: it consists of splitting the length a system into two
blocks (which can initially be of small size), introducing at each step two
additional sites, such that the total system is now dubbed a superblock and
then finding the ground state of this superblock with a diagonalization method
such as the Lanczos algorithm. The resulting set of eigenstates is truncated
so that only the states corresponding to the largest eigenvalues of the
reduced density matrix for a superblock are maintained. These are then used to
build the basis of states corresponding to the new superblocks, which can now
be augmented by two additional sites. This process, called infinite-system
DMRG, is repeated until the desired system size $L$ is reached. It is often
viewed as a warm up cycle for the finite-system DMRG stage, which starts with
the basis found in the infinite-system DMRG stage and continues the growth of
the two blocks in the system, but now keeping the total size fixed such that
each block grows at the expense of the other. The process of each of the
blocks being increased to maximum size as the other shrinks to its minimum
size is called a sweep, and usually several sweeps are required before
convergence in the energy of the system is achieved. The resulting states can
then be utilized to calculate properties such as energy, correlation and order
parameters. While the main focus of DMRG is obtaining the ground state
observables, it can also produce accurate results for the wave function and
for low-lying excited states, within the approximations of the method.

In recent years, DMRG has seen a number of updates and adaptations, including
formulations in terms of Matrix Product States (MPS) (McCulloch, 2007;
Verstraete et al., 2008; Schollwöck, 2011; Wall and Carr, 2012) and the
development of related methods capable of treating dynamics, such as the time-
evolving block decimation (TEBD) algorithm (Vidal, 2003, 2004) and the
adaptive time-depedendent DMRG (Daley et al., 2004; Schollwöck and White,
2006). These tools have been widely employed in the study of discrete many-
body quantum systems, specially in the treatment of Heisenberg-type spin
chains (White and Huse, 1993; Hallberg, 1995; Gobert et al., 2005; Führinger
et al., 2008)

$H=J\sum_{\langle ij\rangle}\mathbf{S}_{i}\cdot\mathbf{S}_{j},$ (92)

and Hubbard-like models (Daul and Noack, 1996; Kühner and Monien, 1998;
Kollath et al., 2005, 2007; Ejima et al., 2011)

$H=-t\sum_{\langle
ij\rangle,\sigma}(c_{i,\sigma}^{\dagger}c_{j,\sigma}+\text{H.c.})+U\sum_{i}n_{i\uparrow}n_{i,\downarrow},$
(93)

and has been useful in the understanding of important phenomena both in the
theory of condensed matter and cold atoms. While originally devised for
studying the low-energy properties of 1D systems, DMRG has also been used to
approach problems in higher dimensions (Nishino, 1995; White, 1996; White and
Scalapino, 1998; Stoudenmire and White, 2012; Orús, 2014), although this is
considerably more difficult due to the scaling with system size beyond 1D.
More recently, a number of finite temperature DMRG studies have also been
performed (Verstraete et al., 2004; Zwolak and Vidal, 2004; Feiguin and White,
2005; White, 2009; Karrasch et al., 2012, 2013), as well as generalizations of
DMRG to continuum models (Verstraete and Cirac, 2010; De las Cuevas et al.,
2018). For further information and for hands-on implementations, we refer to
the different open-source DMRG-based libraries that are currently available
(Bauer et al., 2011; Jaschke and Carr, 2018; Fishman et al., 2020).

###  9.4 The Multi-Layer Multi-Configuration Time-Dependent Hartree Method for
Atomic Mixtures

Another way to approach the quantum many-body problem is to use wave function
propagation methods that address the exponential scaling of the Hilbert space
with the particle number by expanding the many-body state in terms of a
dynamically optimized truncated basis. The Multi-Configuration Time-Dependent
Hartree (MCTDH) family of methods (Meyer et al., 1990; Beck et al., 2000; Lode
et al., 2020) is such an approach for the ab-initio treatment of multi-
dimensional time-dependent problems. MCTDH was initially designed for handling
distinguishable particles with a focus on the multimode vibrational dynamics
of molecules. Early works, dating back about two decades, used it to describe
few-body settings of distinguishable particles, e.g. see (Zöllner et al.,
2006b, a, 2007). To overcome this limitation, MCTDH was consequently extended
to treat bosonic and fermionic systems through the development of MCTDHB (Alon
et al., 2008; Streltsov et al., 2007) and MCTDHF (Zanghellini et al., 2003;
Caillat et al., 2005; Fasshauer and Lode, 2016; Nest, 2006), respectively.
These methods have been successfully applied to study the stationary
properties and, importantly, the correlated quantum dynamics of various
ultracold bosonic and fermionic systems. Notable examples include tunneling
processes (Sakmann et al., 2009; Lode et al., 2012; Schurer et al., 2016),
collective excitations (Schmitz et al., 2013), quench-induced dynamics
(Mistakidis et al., 2015b, 2014; Neuhaus-Steinmetz et al., 2017), periodic
driving of lattice-trapped atoms (Mistakidis et al., 2015a; Mistakidis and
Schmelcher, 2017), the effect of long-range e.g. dipolar interactions
(Chatterjee et al., 2020, 2019; Bera et al., 2020), quantum non-linear
excitations (Streltsov et al., 2011; Krönke and Schmelcher, 2015; Katsimiga et
al., 2018, 2017b; Beinke et al., 2018) and atom-ion systems (Schurer et al.,
2015). Recent developments have aimed tackling the computational complexity of
large systems due to the exponentially growing Hilbert space and exploit
configuration selection schemes for MCTDH (Worth, 2000; Wodraszka and
Carrington, 2016; Wodraszka and Carrington Jr, 2017), a static selection of
the most relevant Hartree products/Slater determinants in MCTDHB(F) (Miyagi
and Madsen, 2013; Lévêque and Madsen, 2017), and dynamical self-adapting
detection of the important number states (Köhler et al., 2019) or the wave
function coefficients (Larsson and Tannor, 2017).

Figure 44: Schematic representation of a three-layer tree diagram containing
(from top to bottom) the top layer, the species layer and the particle layer.
The system corresponds to a three species bosonic mixture with species $A$,
$B$, $C$. The number of species (single-particle) functions are $M_{A}$,
$M_{B}$ and $M_{C}$ ($m_{A}$, $m_{B}$ and $m_{C}$). The primitive degrees-of-
freedom $x_{A}$, $x_{B}$ and ($x_{C}$, $y_{C}$, $z_{C}$) are marked by the
squares and have dimension $\mathcal{M_{A}}$, $\mathcal{M_{B}}$ and
($\mathcal{M}_{C,x}$, $\mathcal{M}_{C,y}$, $\mathcal{M}_{C,z}$). The "+"
symbol within the species nodes designates the fact that the species are
bosonic. Figure from (Cao et al., 2013).

An important extension of the MCTDH family of methods was the integration of
the multi-layer (ML) scheme into the wave function ansatz (Manthe, 2008;
Vendrell and Meyer, 2011; Wang and Thoss, 2003). It is based on a coarse-
graining cascade, where strongly correlated degrees-of-freedom are grouped
together and treated as subsystems that mutually couple to each other (see
Fig. 44). This feature enables one to adapt the wave function ansatz to system
specific intra- and intercomponent correlations and address mesoscopic Bose-
Bose mixures using ML-MCTDHB (Cao et al., 2013; Krönke et al., 2013) or
arbitrary ultracold atomic mixtures such as Bose-Fermi and Fermi-Fermi
mixtures including spinor components (Koutentakis et al., 2019) and different
spatial dimensions within the ML-MCTDHX (Cao et al., 2017) framework. Recent
applications include, but are not limited to, tunneling of mixtures (Keiler et
al., 2020; Theel et al., 2020), phase-separation processes (Mistakidis et al.,
2018; Erdmann et al., 2019), Bose (Mistakidis et al., 2020a) and Fermi
(Mistakidis et al., 2019d) polaron dynamics, quantum dark-bright solitons
(Katsimiga et al., 2017a), mesoscopic charged molecules (Schurer et al., 2017)
and spinor BECs (Mittal et al., 2020). The main concept is to expand the state
of the atomic mixture utilizing dynamically optimized truncated basis states
for the involved components, which are in turn efficiently represented as
MCTDHB(F) wave functions. ML-MCTDHX (Cao et al., 2017) is a flexible
variational method for solving the time-dependent many-body Schrödinger
equation of atomic mixtures while capturing the emergent correlation effects
using a computationally feasible basis size.

In general, an $S$-component ultracold atomic mixture is a composite system
residing in the Hilbert space
$\mathcal{H}^{AB,\dots,S}=\mathcal{H}^{A}\otimes\mathcal{H}^{B}\otimes\dots\otimes\mathcal{H}^{S}$,
where $\mathcal{H}^{\sigma}$ denotes the Hilbert space of the
$\sigma=A,B,\dots,S$ component. To account for the inter-component
correlations of the mixture the many-body wave function can be expanded in
terms of a set of $k=1,2,\dots,D_{\sigma}$ distinct species functions
$\Psi^{\sigma}_{k}(\vec{x}^{\sigma};t)$ for each $\sigma$-component, see also
Fig. 44, as

$\begin{split}\Psi(\vec{x}^{A},\dots,\vec{x}^{S};t)=\sum_{i_{A},i_{B},\dots,i_{S}=1}^{D_{A},D_{B},\dots,D_{S}}A_{i_{A},i_{B},\dots,i_{S}}(t)\Psi^{A}_{i_{A}}(\vec{x}^{A};t)\otimes\Psi^{B}_{i_{B}}(\vec{x}^{B};t)\otimes\dots\otimes\Psi^{S}_{i_{S}}(\vec{x}^{S};t).\end{split}$
(94)

In this expression the
$\vec{x}^{\sigma}=(x_{1}^{\sigma},x_{2}^{\sigma},\dots,x_{N_{\sigma}}^{\sigma})$
are the spatial coordinates and the $A_{i_{A},i_{B},\dots,i_{S}}(t)$ are the
time-dependent expansion coefficients offering a measure of the inter-
component correlations. Also, the functions $\\{\Psi_{k}^{\sigma}\\}$ form an
orthonormal set of $N_{\sigma}$-body wave functions in a subspace of
$\mathcal{H}^{\sigma}$.

For $S>2$ the numbers $D_{\sigma}$ of the species functions do not have to
necessarily coincide. However, in the case of a binary mixture considering
$D_{A}\neq D_{B}$ does not improve the accuracy of the wave function ansatz
since according to the Schmidt decomposition for bipartite systems (Horodecki
et al., 2009) the maximal strength of inter-component correlations depends
solely on $\min\\{D_{A},D_{B}\\}$. Thus, we can set $D_{A}=D_{B}\equiv D$
without loss of generality and the ansatz in Eq. (94) reduces to a truncated
Schmidt decomposition of rank $D$, namely
$\Psi(\vec{x}^{A},\vec{x}^{B};t)=\sum_{k=1}^{D}\sqrt{\lambda_{k}(t)}~{}\Psi^{A}_{k}(\vec{x}^{A};t)\Psi^{B}_{k}(\vec{x}^{B};t)$.
Here, the spectral decomposition
$A_{i_{A}i_{B}}(t)=\sum_{k=1}^{D}U_{i_{A}k}^{-1}(t)\sqrt{\lambda_{k}(t)}U_{ki_{B}}(t)$
of the expansion coefficients $A_{i_{A}i_{B}}$ in Eq. (94) has been used where
$\sqrt{\lambda_{k}(t)}$ are the Schmidt weights. These are referred to as the
natural species populations of the $k$-th species function $\Psi^{\sigma}_{k}$
of the $\sigma$-component since they correspond to the eigenvalues of the
species’ reduced density matrix

$\displaystyle\rho^{N_{\sigma}}(\vec{x}^{\sigma},\vec{x}^{\prime\sigma};t)=\int
dx^{\sigma^{\prime}}_{1}dx^{\sigma^{\prime}}_{2}\dots
dx^{\sigma^{\prime}}_{N_{\sigma^{\prime}}}\Psi^{*}(\vec{x}^{\sigma},\vec{x}^{\sigma^{\prime}};t)\Psi(\vec{x}^{\prime\sigma},\vec{x}^{\sigma^{\prime}};t),$
(95)

where $\sigma\neq\sigma^{\prime}$. The natural species populations allow one
to quantify the presence of inter-component correlations and entanglement, and
if multiple eigenvalues of $\rho^{N_{\sigma}}$ are macroscopically occupied,
i.e. ‘on the order of the system size’, then the system is considered species
entangled, whereas otherwise it is not entangled. When at least two distinct
$\lambda_{k}(t)$’s are finite the many-body state cannot be expressed as a
direct product of two states stemming from $\mathcal{H}^{A}$ and
$\mathcal{H}^{B}$ and thus the system is entangled (Roncaglia et al., 2014).
In this sense, the deviation $1-\lambda_{1}(t)$ quantifies the degree of
system entanglement.

Next, in order to include intra-component correlations into the many-body wave
function ansatz each species function $\Psi^{\sigma}_{k}(\vec{x}^{\sigma};t)$
is expressed as a linear superposition of time-dependent number-states
$|\vec{n}(t)\rangle^{\sigma}$ with time-dependent expansion coefficients
$B^{\sigma}_{k;\vec{n}}(t)$:

$|\Psi_{k}^{\sigma}(t)\rangle=\sum_{\vec{n}}B^{\sigma}_{k;\vec{n}}(t)|\vec{n}(t)\rangle^{\sigma}~{}~{}{\rm
for}~{}~{}k=1,2,\dots,D_{\sigma}.$ (96)

Each number state $|\vec{n}(t)\rangle^{\sigma}$ corresponds to a permanent for
bosons and a determinant for fermions. It is constructed upon $d_{\sigma}$
time-dependent variationally optimized single-particle functions (SPFs)
$\left|\phi_{z}^{\sigma}(t)\right\rangle$, with $z=1,2,\dots,d_{\sigma}$ and
occupation numbers $\vec{n}=(n_{1},\dots,n_{d_{\sigma}})$. The occupations
$n_{i}$ are integers, and in particular $n_{i}\in\\{0,1\\}$ for fermions, and
determine the number of particles in a certain SPF while satisfying the
restriction $\sum_{i=1}^{d_{\sigma}}n_{i}=N_{\sigma}$. This implies that
$\Lambda_{\sigma}=\binom{N_{\sigma}+d_{\sigma}-1}{d_{\sigma}-1}$
[$\Lambda_{\sigma}=\binom{d_{\sigma}}{N_{\sigma}}$] different number states
exist for $N_{\sigma}$ bosons (fermions). Furthermore, the SPFs are expanded
with respect to a time-independent primitive basis of dimension $M_{pr}$,
which could, for example, be a discrete variable representation. The
characteristic structure of the wave function expansion at the different
layers for a three component bosonic mixture is provided in Fig. 44. We note
that these multi-layer concepts can be extended to more than one spatial
dimension, while still handling multiple internal components (Bolsinger et
al., 2017).

In order to address the time-evolution of the ($N_{A}+N_{B}+\dots+N_{S}$)-body
wavefunction $\Psi(\vec{x}^{A},\vec{x}^{B}\dots,\vec{x}^{S};t)$ obeying the
Schrödinger equation of interest under the influence of the Hamiltonian $H$,
one can numerically determine the ML-MCTDHX equations of motion (Cao et al.,
2017) based upon a variational principle such as the Lagrangian (Broeckhove et
al., 1988), McLachlan (McLachlan, 1964) or the Dirac-Frenkel (Dirac, 1930;
Frenkel, 1934) principle for the generalized ansatz of Eqs. (94) and (96). As
a case example, referring to a binary Bose-Fermi mixture of species B and F,
these equations consist of a set of $D^{2}$ ordinary (linear) differential
equations of motion for the coefficients $\lambda_{k}(t)$, coupled to a set of
$D[\binom{N_{B}+d_{B}-1}{d_{B}-1}+\binom{d_{F}}{N_{F}}]$ non-linear
integrodifferential equations for the species functions, and $d_{B}+d_{F}$
non-linear integrodifferential equations for the SPFs (Cao et al., 2017). The
main control parameters of the ML-MCTDHX are $D_{\sigma}$ and $d_{\sigma}$,
and it is possible to use the scheme in different orders of approximation. As
such when choosing $d_{\sigma}=M_{pr}$ and
$D_{\sigma}=\binom{N_{\sigma}+d_{\sigma}-1}{d_{\sigma}-1}$ for bosons
[$D_{\sigma}=\binom{d_{\sigma}}{N_{\sigma}}$ for fermions], the wave function
ansatz is of full configuration-interaction type thus covering the full
Hilbert space of the system. In the other extreme case of neglecting all
correlations by setting $D_{\sigma}=d_{\sigma}=1$ for bosons
($d_{\sigma}=N_{\sigma}$ for fermions), the ML-MCTDHX equations of motion
reduce to the set of coupled Gross-Pitaevskii (Hartree-Fock) equations of
motion (Pitaevskii and Stringari, 2016). For concrete cases with systems that
do not feature very strong correlations, it is typically sufficient to
consider a limited number of species functions and orbitals, and hence to
obtain a very efficient representation of the many-body state. Thus,
$D_{\sigma}$ and $d_{\sigma}$ need to be chosen according to the dominant
correlations of the system under consideration such that the many-body
expansion becomes efficient, i.e. it involves as few coefficients as possible,
while properly capturing the correlations.

A major challenge for many-body simulations is to achieve a desirable degree
of convergence. Within the ML-MCTDHX framework one has to carefully inspect
the truncation order of the total Hilbert space which is determined by the
numerical configuration space $C=(\\{D_{\sigma}\\};\\{d_{\sigma}\\};M_{pr})$.
This is done by systematically increasing $D_{\sigma}$, $d_{\sigma}$, and
$M_{pr}$, and comparing the results for the observables of interest.
Additionally, in limited cases where there are analytical estimations for
certain observables, e.g. the center-of-mass variance (Klaiman and Alon, 2015;
Alon, 2019), they can be employed to test the convergence of the ML-MCTDHX
simulations (Katsimiga et al., 2017a). Finally, we note that the predictions
of ML-MCTDHX have been benchmarked against various numerical data, e.g., for
two particles in a harmonic oscillator (Gwak et al., 2021), for the
generalized Gross-Pitaevskii framework with the Lee-Huang-Yang contribution
Mistakidis et al. (2021), for systems with impurities Koutentakis et al.
(2022); Brauneis et al. (2022), and for spin-chain models (Koutentakis et al.,
2019). In general, ML-MCTDHX behaves as most other ab-initio methods and the
degree of convergence in general depends on the observable of interest and the
system under consideration. Indeed, lower order quantities such as the energy
of the system or densities display a relatively faster convergence as compared
to higher-order observables like the $N_{\sigma}$-body reduced density matrix.

###  9.5 Spin-Chain methods

As discussed in Section 4.7, strongly interacting particles with two or more
internal degrees of freedom can be mapped onto a spin-chain Hamiltonian whose
coupling coefficients are determined by the external trapping geometry. This
subsection illustrates the mapping for a two-component system following
(Volosniev et al., 2015b), see also (Deuretzbacher et al., 2014). For the sake
of discussion, it is assumed that $g_{A}=g_{B}$, allowing for the following
form of the Hamiltonian from Eq. (28)

$\displaystyle H=$
$\displaystyle\sum_{n=1}^{N_{A}}H_{0}(x_{n})+\sum_{n=1}^{N_{B}}H_{0}(y_{n})$
$\displaystyle+\kappa g\sum_{n<m}^{N_{A}}\delta(x_{n}-x_{m})+\kappa
g\sum_{n<m}^{N_{B}}\delta(y_{n}-y_{m})+g\sum_{n=1}^{N_{A}}\sum_{m=1}^{N_{B}}\delta(x_{n}-y_{m}),$
(97)

where $H_{0}(x)=-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial x^{2}}+V(x)$
is the single particle Hamiltonian. In the impenetrable limit ($1/g=0$), the
eigenstates of the system can be written as

$\Psi=\sum_{k=1}^{L(N_{A},N_{B})}a_{k}P_{k}\Phi_{F}(\\{x_{i},y_{j}\\}),$ (98)

where $\Phi_{F}(\\{x_{i},y_{j}\\})$ describes the state of $N=N_{A}+N_{B}$
spinless fermions for the ordering of particles
$x_{1}<...<x_{N_{A}}<y_{1}<...<y_{N_{B}}$. $P_{k}$ is a permutation operator
acting on the coordinates. There are

$L(N_{A},N_{B})=\frac{(N_{A}+N_{B})!}{N_{A}!N_{B}!}$ (99)

combinations, and thus independent coefficients $a_{k}$, which lead to
$L(N_{A},N_{B})$ degenerate states at $1/g=0$.

The states in the vicinity of $1/g=0$ should extremize the slope of the energy
near the fermionization limit, which is obtained from the Hellmann-Feynman
theorem,

$\displaystyle\frac{\partial E}{\partial
g}=\frac{\sum_{n=1}^{N_{A}}\sum_{m=1}^{N_{B}}\langle\delta(x_{n}-y_{m})\rangle+\kappa\sum_{n<m}^{N_{A}}\langle\delta(x_{n}-x_{m})\rangle+\kappa\sum_{n<m}^{N_{B}}\langle\delta(y_{n}-y_{m})\rangle}{\langle\Psi|\Psi\rangle},$
(100)

where $\langle\delta(x)\rangle\equiv\langle\Psi|\delta(x)|\Psi\rangle$. The
numerator in this expression can be calculated using the boundary conditions
due to the zero-range interaction. For particles of the same type, the
conditions are given as

$\displaystyle\left(\frac{\partial\Psi}{\partial
x_{i}}-\frac{\partial\Psi}{\partial
x_{i^{\prime}}}\right)\Bigg{\rvert}_{x_{i}-x_{i^{\prime}}=0^{-}}^{x_{i}-x_{i^{\prime}}=0^{+}}=2\frac{\kappa
gm}{\hbar^{2}}\Psi(x_{i}=x_{i^{\prime}}),$
$\displaystyle\left(\frac{\partial\Psi}{\partial
y_{i}}-\frac{\partial\Psi}{\partial
y_{i^{\prime}}}\right)\Bigg{\rvert}_{y_{i}-y_{i^{\prime}}=0^{-}}^{y_{i}-y_{i^{\prime}}=0^{+}}=2\frac{\kappa
gm}{\hbar^{2}}\Psi(y_{i}=y_{i^{\prime}}),$ (101)

and for the interaction between $A$ and $B$ particles one has

$\left(\frac{\partial\Psi}{\partial x_{i}}-\frac{\partial\Psi}{\partial
y_{j}}\right)\Bigg{\rvert}_{x_{i}-y_{j}=0^{-}}^{x_{i}-y_{j}=0^{+}}=2\frac{gm}{\hbar^{2}}\Psi(x_{i}=y_{j}).$
(102)

Equations (100), (101) and (102) with the wave function proposed in Eq. (98)
lead to the expression for the energy in the vicinity of $1/g=0$

$\displaystyle
E=E_{F}-\frac{1}{g}\frac{\sum_{i=1}^{N-1}\alpha_{i}\mathcal{A}_{i}}{\sum_{k=1}^{L(N_{A},N_{B})}a_{k}^{2}},$
(103)

where $E_{F}$ is the eigenenergy of a system of $N_{A}+N_{B}$ spinless
fermions, and

$\mathcal{A}_{i}=\sum_{k=1}^{L(N_{A}-1,N_{B}-1)}(a_{i|k}^{AB}-a_{i|k}^{BA})^{2}+\frac{2}{\kappa}\sum_{k=1}^{L(N_{A}-2,N_{B})}(a_{i|k}^{AA})^{2}+\frac{2}{\kappa}\sum_{k=1}^{L(N_{A},N_{B}-2)}(a_{i|k}^{BB})^{2}.$
(104)

The coefficient $a_{i|k}^{ss^{\prime}}$ corresponds to $a_{k}$ of Eq. (98)
that multiplies $\Phi_{F}$ with $s$ particle at the position $i$ and
$s^{\prime}$ particle at the position $i+1$. There is no need to specify the
order of other particles, in fact the sum over $k$ is the sum over the
ordering of all other particles. The coefficients $\alpha_{i}$ depend only on
the shape of the trapping potential and on the total number of particles
$N=N_{A}+N_{B}$. In particular, they are the same for bosonic and fermionic
systems and are given by

$\alpha_{i}=\frac{\hbar^{4}}{m^{2}}\frac{\int_{x_{1}<x_{2}\cdots<x_{N-1}}dx_{1}\cdots
dx_{N-1}\Big{|}\frac{\partial\Phi_{F}(x_{1},\cdots,x_{i},\cdots,x_{N})}{\partial
x_{N}}\Big{|}^{2}_{x_{N}=x_{i}}}{\int_{x_{1}<x_{2}\cdots<x_{N-1}}dx_{1}\cdots
dx_{N}|\Phi_{F}(x_{1},\cdots,x_{i},\cdots,x_{N})|^{2}}.$ (105)

The mapping of a strongly interacting system onto a spin chain is established
upon observation that the energy functional in Eq. (103) can be obtained with
the wave-function ansatz

$\Psi_{s}=\sum_{k=1}^{L(N_{A},N_{B})}a_{k}P_{k}|\uparrow_{1}...\uparrow_{N_{A}}\downarrow_{1}...\downarrow_{N_{B}}\rangle$
(106)

for the spin-chain Hamiltonian

$H_{s}=\left(E_{F}-\frac{\kappa+1}{g\kappa}\sum_{i=1}^{N_{A}+N_{B}-1}\alpha_{i}\right)\mathbf{I}+\frac{1}{g}\sum_{i=1}^{N_{A}+N_{B}-1}\alpha_{i}\left[\mathbf{P}_{i,i+1}-\frac{1}{\kappa}\sigma_{z}^{i}\sigma_{z}^{i+1}\right],$
(107)

where $\mathbf{P}_{i,i+1}$ is the operator that swaps the spins at the $i$th
and ($i+1$)th sites, $\sigma_{z}^{i}$ is the third Pauli matrix that acts on
the site $i$ and $\mathbf{I}$ is the identity operator. This correspondence
implies that to find the set $\\{a_{k}\\}$ that minimizes the energy in Eq.
(103), it is enough to find the eigenstates of the Hamiltonian $H_{s}$.

The non-trivial part of applying the mapping is the computation of the
$\alpha_{i}$. Given the external potential $V$, it is typically easy to
compute a set of $N=N_{A}+N_{B}$ one-body eigenstates and then construct the
Slater determinant, $\Phi_{F}$, from these $N$ states. However, standard
techniques for multidimensional integration allow one to calculate the
multiple integral in Eq. (105) only for $N\lesssim 15$ (Volosniev et al.,
2014a; Levinsen et al., 2015; Yang et al., 2015; Yang and Pu, 2016). It is
hard, if not impossible, to extend multidimensional integration to a large
number of particles ($N\simeq 15$), due to the curse of dimensionality, see,
e.g., (Hinrichs et al., 2014) and references therein.

A method to tackle systems with $N>15$ is highly desirable to study a number
of research questions, such as few- to many-body transition, and finite-size
effects in mesoscopic cold-atom experiments. Luckily, Eq. (105) can be
rewritten as a sum of one-dimensional integrals (Loft et al., 2016a), which
allows one to calculate the $\alpha_{i}$ for systems with $N\lesssim 35-60$
(Loft et al., 2016a; Deuretzbacher et al., 2016). An open source code CONAN
(Coefficients of One-dimensional $N$-Atom Networks) for such calculations is
also available (Loft et al., 2016a).

To go to even larger systems (e.g., with $N\simeq 60$), one could use the
existing approximate expressions of $\alpha_{i}$ for different confining
potentials. For large systems in external traps that change weakly on the
length scale given by the density of spinless fermions, the local density
approximation (Deuretzbacher et al., 2014) can provide one with an accurate
expression

$\alpha_{i}=\frac{\hbar^{4}\pi^{2}n_{TF}^{3}(z_{i})}{3m^{2}},$ (108)

where $n_{TF}$ is the Thomas-Fermi density profile, and $z_{i}$ is the center
of mass of the $i$th and the ($i+1$)th particle. The expression in Eq. (108)
can be anticipated from the results obtained for homogeneous quantum wires
(Matveev, 2004). Its accuracy has been checked for various trapping geometries
(Deuretzbacher et al., 2014; Marchukov et al., 2016; Yang and Pu, 2016), where
it was shown that Eq. (108) is qualitatively correct even for a handful of
particles, provided that the trap changes weakly on the length scale given by
the density of the gas. Although, quantitative differences between the exact
values and the ones obtained with Eq. (108) might lead to noticeable effects
in time dynamics for a small numbers of particles (Marchukov et al., 2016),
Eq. (108) is a powerful starting point for understanding the relationship
between the trapping potential and the coefficients $\alpha_{k}$. In reverse-
engineering applications, Eq. (108) provides an initial guess for trapping
potentials that may lead to perfect state transfer (Loft et al., 2016b) or
that may realize desired entanglement Hamiltonians (Barfknecht et al., 2021b).

For a harmonic trap, the values of the $\alpha_{i}$ have been approximated by
(Levinsen et al., 2015)

$\alpha_{i}\simeq C_{N}i(N-i),$ (109)

where the coefficient $C_{N}$ depends only on the total number of particles.
For a given system with $N$ fixed, it sets the overall energy scale; its value
is not important for calculations of particle-particle correlations. The
values of $C_{N}$ including its asymptotic expansion are presented in
(Massignan et al., 2015). Equation (109) is not exact, but it is surprisingly
accurate (Levinsen et al., 2015; Yang et al., 2015; Laird et al., 2017). Its
simplicity might imply some hidden symmetries for strongly interacting systems
in harmonic traps.

For a box potential262626(Duncan et al., 2017) observed similar physics in an
optical lattice potential for certain fillings., one expects a Bethe-ansatz
solvable spin chain. This expectation is supported by indications that the
exchange coefficients do not depend on the position. This was seen numerically
for the ground (Marchukov et al., 2016) and excited states (Barfknecht et al.,
2021a). For the ground state, this observation can be supported by the Bethe
ansatz results (Pan et al., 2017). Assuming that the coefficients are indeed
identical, one can show that272727Note that this expression could be used to
derive Eq. (108). $\alpha_{i}=2\hbar^{2}E_{F}/(mL)$ (Volosniev, 2017;
Barfknecht et al., 2021a). This simple form of the coefficients allows for the
investigation of static and dynamic properties of strongly interacting systems
in a box using homogeneous spin chains (Barfknecht et al., 2021a). In
addition, it provides a starting point for analyzing systems with barriers, as
discussed in (Iversen et al., 2020).

##  10 Beyond contact interactions

(Almost) all physics discussed in the previous chapters is based on the
assumption that two-body interactions can be modelled using delta-function (or
contact) potentials, $g\delta(x_{1}-x_{2})$ – the standard approach for cold-
atom systems in 1D, for a review see (Yurovsky et al., 2008). Unlike higher
dimensions, see (Braaten and Hammer, 2006), the zero-range two-body
interaction leads to a well-defined many-body problem in 1D. Therefore, from a
purely mathematical point of view, there is no need for theoreticians to
consider other (more fundamental) interactions. However, from the physical
standpoint, the potential, $g\delta(x)$, has its limitations, and cannot be
used to describe a number of important phenomena. For example, delta-function
interactions cannot bind spinless fermions (due to the Pauli principle) or
lead to Borromean and Brunnian systems282828Borromean (Brunnian) systems are
bound states of more than two particles in which two-body (all) subsystems are
unbound, see, e.g., (Baas et al., 2014). In 3D, many properties of such states
can be explored using zero-range potentials assuming a three-body parameter
(Jensen et al., 2004; Braaten and Hammer, 2006). In 1D and 2D, finite-range
potentials are necessary for Borromean and Brunninan systems since any
attractive potential supports a two-body bound state (Landau and Lifshitz,
1977).. These observations imply that delta-function potentials cannot be used
to model physics of, for instance, unconventional superconductors or
superfluids [such as classic helium-3 systems (Leggett, 1972)], or of exotic
systems in which few-body settings play the role of building blocks.

A modification of the delta-function potential is therefore necessary for
studying certain types of new physics. This necessity motivates theoretical
and experimental studies beyond contact interactions. There are a few well-
known systems with such interactions, and this section briefly reviews some of
them. Our aim here is to provide a succinct description of some rapidly
developing topics lying beyond the main focus of the present review.

###  10.1 Finite-range corrections in atom-atom interactions

In cold-atom physics, atom-atom interactions are usually considered to be of
short range character. The validity of this approximation has been established
for many systems by analyzing various experimental data with theoretical
models based upon $g\delta(x)$, where $g$ then describes only some integral
information about the actual potential, see, e.g., Fig. 41. It turns out,
however, that finite-range corrections are important for narrow Fano-Feshbach
resonances, which are closed-channel dominated (Friedrich and Friedrich, 2006;
Greene et al., 2017) yielding an energy-dependent scattering length (Chin et
al., 2010).

A similar scenario may arise for collisions in the presence of quasi-1D set-
ups, where the transverse spatial directions are tightly confined as compared
to the axial one (i.e., $\omega_{x},\omega_{y}\gg\omega_{z}$). Then, it might
be beneficial to incorporate finite-range corrections into the coupling
constant $g$ regardless of the type of Fano-Feshbach resonances that are
employed to tune the two-body atomic interactions. To model two-body physics
of such systems, one can for example use the energy-dependent scattering
length

$a(k)\simeq\frac{a}{1-\frac{1}{2}k^{2}ar_{e}},$ (110)

where $k$ is the two-body momentum in the relative two-body frame given by the
relation $k=\sqrt{2\mu E/\hbar^{2}}$ with $\mu$ indicating the reduced mass of
the particles and $E$ their total collisional energy. Also, $r_{e}$ is the
effective range correction and $a$ denotes the energy-independent scattering
length parameter.

In quasi-1D trapping potentials the relation (110) was used in (Naidon et al.,
2007; Yurovsky et al., 2008) (see (Blume and Greene, 2002) for relevant 3D
discussions) to construct models that allow for the description of resonant
phenomena, such as confinement-induced resonances (Yurovsky et al., 2008;
Dunjko et al., 2011). Note that many-body problems with energy-dependent
coupling coefficients require a modification of the commonly used numerical
approaches, and therefore, their physics is less explored.

###  10.2 Odd-wave interactions

Simulations of ‘unconventional’ many-body physics in 1D may be based upon
three-dimensional $p$-wave physics, which leads to interactions between two
spin-polarized fermions. The strength of the interactions can be tuned using
external fields, see, e.g., (Granger and Blume, 2004; Peng et al., 2014).

Severe atom losses complicate studies close to $p$-wave resonances in 3D, see,
e.g., (Regal et al., 2003; Luciuk et al., 2016). However, recent theoretical
and experimental work (Günter et al., 2005; Zhou and Cui, 2017; Kurlov and
Shlyapnikov, 2017; Chang et al., 2020) may provide a background for future
experimental realization of strongly interacting one-dimensional systems with
odd-wave interactions. As an example, Figure 45 presents experimental results
for the dependence of the peak three-body loss coefficient, $L_{3}$, with
respect to the lattice depth, $V_{L}$, in 6Li (Chang et al., 2020). Only a
minor dependence on the confinement strength is visible, and (Chang et al.,
2020) suggest that suppressing loss in heavier fermions such as 40K, 161Dy and
167Er might be even more promising. Those systems may then allow one to study
few- and many-body statics and dynamics with odd-wave interactions.

Figure 45: Peak three-body loss coefficient for fermions interacting via a
$p$-wave resonance as a function of of the lattice depth. Figure from
supplemental material of (Chang et al., 2020).

To describe systems with odd-wave interactions, one requires a beyond-contact
potential, since the term $g\delta(x)$ acts only in the even-wave channel due
to the symmetry $\delta(-x)=\delta(x)$. A possible route to describe odd-wave
interactions is via the pseudopotential

$v(x)=g_{P}\delta^{\prime}(x)\hat{O}_{\pm},\qquad\hat{O}_{\pm}f(x)=\frac{1}{2}\left(f^{\prime}(0+)+f^{\prime}(0-)\right),$
(111)

which defines the odd-wave boundary conditions for the wave function
(Girardeau and Olshanii, 2004; Girardeau et al., 2004). Note that this zero-
range approximation to odd-wave potentials in 1D requires regularization
(Valiente and Zinner, 2015; Cui, 2016), which is similar to zero-range
$s$-wave problems in 3D ((Braaten and Hammer, 2006)). A theoretical approach
that includes beyond-zero-range physics and does not require a regularization
can be based upon the finite-range potential (Veksler et al., 2014):

$v(x)=g_{0}\delta(x)+g_{1}[\delta(x-d)+\delta(x+d)],$ (112)

where $g_{0},g_{1}$ and $d$ are parameters that determine the strength and the
range of the interactions.

The experimental advances on odd-wave interactions have motivated
theoreticians to study stability and phase diagrams of many-body systems with
odd-wave interactions (Hao et al., 2007; Imambekov et al., 2010; Pan et al.,
2018). Theoretical tools for such studies are the Bethe ansatz, which can be
used under certain conditions, in particular when the Bose-Fermi mapping is
applicable (Cheon and Shigehara, 1999; Valiente, 2020, 2021). To tackle many-
body problems that cannot be solved by the Bethe ansatz, one can rely on
numerical methods discussed in the previous section as well as on variational
approaches, see, e.g., (Kościk and Sowiński, 2020). For strongly interacting
systems one can employ effective spin-chain formulations (Hu et al., 2016b),
and it is possible to use a mapping of a continuous model onto a lattice
representation (Muth et al., 2010). Several such approaches inspired by
effective field theory (Hammer and Furnstahl, 2000) and few-body routes to the
physics of liquids have been introduced recently (Valiente and Öhberg, 2016;
Valiente, 2020; Morera et al., 2021b). The physics of the beyond Lieb-Liniger
model that corresponds to the interaction in Eq. (112) is discussed in
(Veksler and Fishman, 2016; Jachymski et al., 2017).

Generalizations of the $s$-wave confinement induced resonances [see also
Section 1.4] to higher partial wave two-body interactions is considered in
(Giannakeas et al., 2011, 2012; Hess et al., 2015). Specifically, it is
demonstrated that the inclusion of higher orbital angular momentum, i.e.,
$\ell>0$, leads to a series of confinement-induced resonances coupled to each
other due to the trapping potential. The energy dependence of such resonances
enjoys several universal attributes, such that for deep enough interatomic
potentials only four $\ell$-wave confinement-induced resonances are important,
i.e., two for identical bosons and two for identical fermions (Hess et al.,
2014). The $\ell$-wave confinement-induced resonances can be important in
dipolar gases (Giannakeas et al., 2013) or in ultracold gases in the vicinity
of high-$\ell$ Fano-Feshbach resonances. Indeed, recent experimental advances
reported the existence of high-partial wave Fano-Feshbach resonances where
open and closed channels are dominated by $\ell\geq 2$ angular momentum in
binary mixtures of 87Rb–85Rb atoms (Cui et al., 2017). Furthermore, the
experimental observation of the first Bose gas near a $d$-wave shape resonance
introduced alternative possibilities where $\ell$-wave confinement-induced
resonances can play a pivotal role (Yao et al., 2019).

###  10.3 Dipole-dipole interactions

A conceptually different approach to simulate beyond-zero-range physics is
based on long-range interactions292929In the context of the present review,
long range means that the potential has a tail which decays slower than the
tail of the van der Waals potential, which characterizes the long-range
behavior of the interaction among neutral atoms.. Note that the one-
dimensional geometry allows one to solve analytically problems with certain
long-range potentials, notably with $1/x^{2}$ potentials (Calogero, 1971).
Those potentials can describe cold-atom experiments only in a relatively crude
approximation, and will not be discussed here. Instead, the focus will be on
the experimentally more relevant case of dipole-dipole potentials that occur
naturally in systems of cold molecules or atoms with large magnetic moment
(e.g., Cr, Dy, Er) (Lahaye et al., 2009; Baranov et al., 2012).

The systems discussed in this subsection can be used to engineer long-range
physics, which is often necessary for cold-atom simulations. For example,
magnetic atoms (Kotochigova, 2014; Norcia and Ferlaino, 2021) were used to
produce BEC of dipoles (Griesmaier et al., 2005; Lu et al., 2011; Aikawa et
al., 2012; Petter et al., 2019) and in particular to realize supersolid states
of matter (Böttcher et al., 2019; Sohmen et al., 2021). A supersolid refers to
a peculiar state of matter where a frictionless flow of a superfluid and
crystal periodic density modulations of a solid coexist (Boninsegni and
Prokof’ev, 2012; Leggett, 1970; Yukalov, 2020). The initial experimental
realization of supersolids in ultracold atom systems were reported in BECs
coupled to light fields created by optical cavities (Léonard et al., 2017) or
featuring spin-orbit coupling (Li et al., 2017). We also note that cold
molecules can have a very large dipole moment, and may lead to other exciting
states of matter. However, some further experimental development is needed as
molecular systems typically suffer from losses more significant that atomic
ones, see (Bohn et al., 2017) for a recent review of the topic.

Systems with dipole-dipole interactions does not have Bethe ansatz solutions
in a continuum, and the $1/x^{3}$ divergence at the origin complicates
calculations in much of the standard numerical machinery303030 It is worth
noting that in the context of models that can be solved by the Bethe ansatz
methods, there is a family of integrable models that describe dipolar atoms in
multi-well systems (Tonel et al., 2015; Ymai et al., 2017) for different
physical applications (Tonel et al., 2020; Wittmann Wilsmann et al., 2018;
Grün et al., 2020, 2022). Also, a variational Bethe Ansatz approach for
dipolar one-dimensional bosons, based on analytic expressions for the ground
state energy and structure factor of the Lieb-Liniger model, has been
discussed (De Palo et al., 2020). From the experimental side, the far-from-
equilibrium dynamics of a strongly interacting nearly-integrable system as it
is tuned away from integrability has been investigated in a dipolar quantum
Newton’s cradle consisting of highly magnetic dysprosium atoms (Tang et al.,
2018). Dysprosium atoms have also been used to create strongly correlated
prethermal states by topological pumping to stabilize a strongly attractive
dipolar system in 1D (Kao et al., 2021). . Luckily, in quasi-one-dimensional
geometries a dipole-dipole potential is regularized at short distances (Sinha
and Santos, 2007; Deuretzbacher et al., 2010; Zinner et al., 2011) and can be
written in the form

$V(x)=\frac{D^{2}}{l^{3}}\left[\frac{1}{4}F\left(\frac{x}{l}\right)-\frac{2}{3}\delta\left(\frac{x}{l}\right)\right],$
(113)

where $D$ is the dipole moment, $l$ is the length associated with the
confinement, and the bounded function $F(u)$ has the expression

$F(u)=-2u+\sqrt{2\pi}(1+u^{2})e^{-\frac{u^{2}}{2}}\mathrm{Erfc}(u/\sqrt{2}).$
(114)

The potential $V(x)$ is determined by the delta-function potential at the
origin, and decays as $1/|x|^{3}$ as $|x|\to\infty$.

At the level of few-body systems, the interaction given by $V$ allows one to
study (quasi)-crystallization of a harmonically trapped system, in particular
the transition from a fermionized state to a Wigner-like crystal in the ground
state properties (Deuretzbacher et al., 2010; Zöllner, 2011; Kościk, 2017;
Bera et al., 2019) [see (Cremon et al., 2010) for relevant studies in two
dimensions]. An external lattice allows one to effectively tune this
transition, and study the bottom-up assembly of lattice models with next-to-
nearest-neighbor interactions.

Figure 46: One-body density as function of the dipolar coupling $g_{d}$. For
$g_{d}\leq 0.8$ the system forms an SMI state, whose density is mainly
concentrated in the three central wells (superfluid fraction). A larger
$g_{d}$ results in the KCS phase characterised by a two hump structure of the
density profile in the doubly occupied wells. An even stronger dipolar
strength leads to the formation of DWCS with a density configuration
alternating between single and double occupations. (b) Natural populations
$\lambda_{k}$, $k=1,2,\dots,8$ with respect to $g_{d}$. In the SF the
$\lambda_{1}$ is dominant while an increasing $g_{d}$ gives rise to
fragmentation which becomes substantial beyond the KCS state. In all cases,
$N=5$ bosons are trapped in a five well lattice potential with barrier height
$V_{0}=8$. All quantities are given in terms of the recoil energy. Figure is
from (Chatterjee et al., 2020).

The ground-state phase diagram of a few dipolar bosons confined in both
commensurate (Chatterjee et al., 2013) and incommensurate (Chatterjee et al.,
2013; Chatterjee and Lode, 2018; Chatterjee et al., 2019, 2020) finite size 1D
optical lattices has been explored relying on the MCTDHB method, see Section
9.4. It has been shown that various states – the few-body precursors of the
phases in the thermodynamic limit – emerge for increasing dipolar strength,
see also (Pizzardo et al., 2016) for the creation of cat-states in a double-
well. Their existence is a consequence of the competition between the kinetic,
potential, and the involved interaction energies. These states range from a
pure superfluid (SF) for weak dipole-dipole interactions313131The employed
interaction potential has the form
$V_{d}=g_{c}\delta(x_{i}-x_{j})+g_{d}/(\left|x_{i}-x_{j}\right|^{3}+\epsilon)$
with $\epsilon$ being the short-scale cutoff introduced by the transversal
confinement i.e. $\epsilon=a_{\perp}^{3}$ (Deuretzbacher et al., 2010, 2013;
Sinha and Santos, 2007) regularizing the divergence at small separations
$\left|x_{i}-x_{j}\right|\leq a_{\perp}$. Here, $a_{\perp}$ is the transveral
length scale. Note also that for electric (magnetic) dipoles
$g_{d}=\mu_{m}^{2}/(4\pi\epsilon_{0})$ ($g_{d}=\mu_{m}^{2}\mu_{0}/(4\pi)$)
where $\mu_{m}$ refers to the dipole moment and $\epsilon_{0}$ ($\mu_{0}$) is
the vacuum permittivity (permeability)., $g_{d}$, which turns into a Mott-
insulator coexisting with a superfluid (SMI) for larger couplings as can be
seen in Fig. 46(a). Here, due to the hard-wall boundaries the superfluid
fraction resides in the central wells in order to minimize the kinetic energy.
A further increase of the dipolar interaction results in the so-called kinetic
crystal state (KCS) 323232Note that the KCS phase does not exist for larger
system sizes and only the DWCS persists., where the bosons avoid each other
and the density exhibits a two-hump structure in the doubly occupied wells.
The latter correspond to the central sites due to the kinetic energy. For even
stronger dipole-dipole couplings, where the interaction energy overcomes the
kinetic energy, a density wave ordered structure appears, in which next-to-
the-nearest neighboring wells are doubly occupied and their density feature a
two-hump structure. The resultant phase is called a density-wave crystal state
(DWCS). It is characterized by substantial two-body anti-correlations within
each lattice site. This behavior is general for crystal phases (Chatterjee et
al., 2013; Chatterjee and Lode, 2018; Chatterjee et al., 2019). As a further
illustration of the many-body nature of these systems, Fig. 46(b) shows the
fragmentation parameter $\Delta=\sum_{k}(\lambda_{k}/N)^{2}$ ($N$ is the atom
number), which serves as a mesoscopic order parameter characterizing these
phases. In this expression, $\lambda_{k}$ denotes the $k$-th eigenvalue of the
reduced single-particle density matrix (Sakmann et al., 2008). It is apparent
that the formation of a crystal-like state is associated with a population
$\Delta=1/N$ of the eigenvalues, a result that indicates the strongly
correlated origin of these states, see also (Fischer et al., 2015) for the
dependence of the degree of fragmentation in long-range settings on the
spatial dimension. An important point here is that in contrast to the SF to
Mott transition, which is inherently related to the presence of an external
lattice potential, crystalization is a result of the dominant effect of the
dipolar interactions and occurs irrespectively of the boundary conditions and
the one-body potential (Arkhipov et al., 2005; Astrakharchik et al., 2008;
Zöllner, 2011; Bera et al., 2019). Interestingly, an experimental protocol has
been proposed for the detection of the above-described phases (Chatterjee et
al., 2020), based on the particle number distributions obtained through
single-shot absorption images (Bakr et al., 2009; Sherson et al., 2010).

The multi-band tunneling dynamics of two dipolar bosons in a 1D triple well is
studied in (Chatterjee et al., 2013). In particular, the authors connect
transitions between bands to the two boson energy spectrum, which features
several avoided-crossings between eigenstates of different bands. The bosons
are initially prepared in their ground state in the most left well of the
triple well and are then left to evolve experiencing a fixed dipole-dipole
coupling lying close to the aforementioned resonances. Depending on the value
of the dipolar strength, being below or above the location of the resonance, a
different dynamical response is observed. It ranges from pair tunneling to
multifrequency tunneling events whose period is substantially impacted by the
long-range coupling. The characteristics of the dynamics are monitored via the
population of each well and the pair probability in the course of the
evolution while even the participating multi-band number states can be
identified due to the few-body nature of the system. It is also worthwhile to
note that using the Floquet picture a directed selective tunneling pathway of
($N-1$) atoms was found to be excited for dipolar bosons which are initialized
in the middle well of a high frequency shaken triple well potential (Lu et
al., 2014; Luo et al., 2015). Note that for dipolar bosons in a double-well
potential the dominant tunneling mechanisms remain qualitatively the same as
for contact interactions. They range from Josephson oscillations to
macroscopic quantum self trapping as it has been discussed e.g. in (Haldar and
Alon, 2018; Roy et al., 2022; Mazzarella and Dell’Anna, 2013).

Figure 47: Illustration of different complexes that can form in layered
geometries with long-range dipolar interacting particles. The complexes are
relevant for both 1D tubes and 2D layers. Notice that the dipole moments are
perpendicularly oriented in these examples. Figure from (Volosniev et al.,
2012)

Beyond a single-tube geometry, the long-range dipolar interaction motivates
studies of coupling between tubes. The inter-tube coupling is overall
attractive and supports a two-body bound state (Zinner et al., 2011; Volosniev
et al., 2013), which may lead to the existence of few-body bound states
between tubes, and, hence, a new type of quantum matter where few-body bound
states play the role of building blocks (Wang et al., 2006; Armstrong et al.,
2020). There are many types of clusters of few dipoles. Chains can include
only one particle per layer, but more complicated structures contain more that
a single bound particle from one layer, and can create a bifurcated many-
particle bound state in layered systems. Such structures can be created by
tilting the direction of dipoles (Wunsch et al., 2011; Volosniev et al.,
2013), and probed for instance by detecting scattered light (Mekhov, 2013). An
illustration of various complexes that could be possible with long-range
dipolar interactions is shown in Figure 47. Note that while the figure
illustrates the case of 2D layers and perpendicularly oriented dipole moments,
the same complexes are imaginable with 1D tubular geometries. In particular,
if the orientations of the dipole moments are tilted with respect to the
planes/tubes, the relative contributions of attraction and repulsion can be
changed (Zinner et al., 2011).

Figure 48: The peak density in a gas of $N=5000$ atoms with respect to
$l_{\perp}/a$ and $\epsilon_{dd}$ quantifying the regimes of attractive and
repulsive LHY contribution in 1D. For $\epsilon_{dd}>1$ a smooth crossover
from droplets to solitons is realized. The white solid line marks the region
where the LHY contribution is on-average repulsive. Here
$l_{\perp}=\sqrt{\hbar/(M\omega_{\perp})}$ is the transversal oscillator
length with frequency $\omega_{\perp}$. Figure from (Edler et al., 2017)

###  10.4 Dipolar droplets

As discussed in Section 8, quantum droplets may appear in interacting Bose-
Bose mixtures with short-range interactions due to the presence of quantum
fluctuations. In these settings the liquid phase occurs as a result of the
competition between the repulsive intracomponent interactions and the
attractive intercomponent coupling. In the same manner droplet structures may
emerge in single-component dipolar bosonic systems due to the interplay
between the long-range dipole-dipole attraction and the short-range repulsion
(Böttcher et al., 2020). In fact, an ensemble of confined attractive dipoles
is typically unstable and its stability can be achieved by the counteraction
of strong repulsive contact interaction that prevents the collapse of such
settings. However, in the case of an overall mean-field attraction the system
is stabilized only by the presence of quantum fluctuations. Dipolar quantum
droplets were observed experimentally for the first time using 164Dy (Ferrier-
Barbut et al., 2016a; Kadau et al., 2016) soon after their theoretical
prediction in Bose mixtures (Petrov, 2015). Since then they have been realized
using several different atomic species including 133Er (Chomaz et al., 2016),
164Dy (Schmitt et al., 2016; Ferrier-Barbut et al., 2016b) and 162Dy (Böttcher
et al., 2019). One of the conclusions drawn from these experiments is that
beyond mean-field effects might be more important in dipolar condensates than
in those with short-range interactions.

Figure 49: Representative dipolar droplet collisions. Panel (a) depicts the
propagation of a long-lived droplet dimer while panels (b)-(e) present in
($\delta=0$) and out-of-phase ($\delta=\pi$) dynamics for various velocities
(see legends). The initial velocity is $u_{0}$, $\xi(\epsilon_{dd})$ is the
dipolar healing length and $m$ the atom mass. Figure from (Edmonds et al.,
2020).

It has been demonstrated (Edler et al., 2017) that the momentum dependence of
the dipolar interactions in quasi-1D geometries (Sinha and Santos, 2007) leads
to beyond mean-field corrections that differ significantly from those in non-
dipolar counterparts. In particular, besides the difference in the density
dependence of the quantum fluctuations in the presence and absence of long-
range interactions also their sign can be changed due to the crucial role of
the transversal directions in 1D dipolar gases. For a relevant
characterization of the emergent phases see Fig. 48. An important control
parameter for entering different phases in these gases is the relative ratio
between the dipolar and the $s$-wave interaction strengths, namely
$\epsilon_{dd}=\mu_{0}\mu_{m}^{2}/(6\pi l_{\perp}^{2}g_{1D}$). Here, the
effective 1D interaction strength is related to the 3D one via
$g_{1D}=g_{3D}/(2\pi l_{\perp}^{2})$ with the transversal length scale being
$l_{\perp}=\sqrt{\hbar/(M\omega_{\perp})}$, while $\mu_{0}$ is the vacuum
permeability and $\mu_{m}$ the dipole moment. Specifically, for dipolar bosons
the density dependence of the LHY term is $\sim-n_{1D}$ as long as $n_{1D}a\to
0$, while for increasing $n_{1D}a$ it departs from this linear behavior (Edler
et al., 2017). Here, $n_{1D}$ denotes the gas peak 1D density and $a$
represents the 3D $s$-wave scattering length. Recall that in Bose mixtures the
LHY scales as $\sim n_{1D}^{1/2}$. It should be clarified at this point that
this behavior is in sharp contrast to the 3D scenario333333Similarly, three-
body losses exhibit a distinct density dependence for quasi-1D dipolar and
non-dipolar Bose gases (Edler et al., 2017), a behavior that is not
encountered in 3D. where the density dependence of the corresponding LHY
contribution in both dipolar gases and Bose-Bose mixtures is $\sim
n_{3D}^{3/2}$ (Böttcher et al., 2020; Petrov, 2015). Moreover, it was found in
(Smith et al., 2021; Bisset et al., 2021) that more complex phases forming
miscible and immiscible droplet states can be achieved in binary magnetic
gases, namely two-component dipolar bosonic gases. The excitation spectrum of
self bound dipolar droplets using Bogoliubov-de-Gennes linearization analysis
has been examined in (Böttcher et al., 2020; Baillie et al., 2017; Lima and
Pelster, 2012). A main result was the identification of modifications e.g. in
the frequencies of the collective modes caused by the inclusion of quantum
fluctuations in the presence of dipolar interactions. Furthermore, by
monitoring the behavior of the low-lying excitation modes it was found that
for increasing atom number the droplets transit from being compressible to
incompressible.

The collisional properties of two counterpropagating 1D dipolar droplets with
the constant velocities $\pm v_{0}$ and the phase difference $\delta$ are
explored in (Edmonds et al., 2020). It is shown that the droplets undergo
inelastic collisions, that their post-collision shapes and sizes depend
crucially on $\delta$, and that excitations build up on the droplets in the
form of sound waves. It is also shown that, for a fixed phase difference,
varying the initial velocity gives rise to a rich droplet response.
Characteristic examples are presented in Fig. 49. One can readily show that
two initially out-of-phase droplets produce a long-lived droplet dimer. Also,
by considering different finite initial velocities it is possible to conclude
that in-phase collisions lead to the droplet fission (Dingwall et al., 2018)
associated with the generation of several droplet states after collision [see
Fig. 49 (b), (d)]. Turning to out-of-phase collisions, one can observe the
formation of two droplets [see Fig. 49 (c), (e)] accompanied by the emission
of sound, which is more pronounced for large initial velocities.

###  10.5 Atom-ion interactions

In recent years, there has been a new push to combine cold atoms with trapped
ions in order to leverage the advantages of both platforms, for a review see
(Tomza et al., 2019). In the context of many of the previous sections
concerned with various types of impurities in backgrounds, we could be tempted
to consider the atom-ion setup with the typical picture of a single ion in a
gas of atoms with the ion playing the role of an impurity. However, the long-
range nature of the atom-ion interactions may complicate the usual polaron
picture and requires a more careful treatment (Côté et al., 2002;
Astrakharchik et al., 2021; Christensen et al., 2021).

The atom-ion interaction of a trapped ion immersed in a cold atomic gas falls
off with the inverse distance to the fourth power, although as with the
dipolar interactions discussed above, the short-range part of the potential
needs regularization (see, e.g., (Krych and Idziaszek, 2015)). Once this is
taken into account, there are many opportunities to study novel physics in the
atom-ion setups, including collision dynamics (Idziaszek et al., 2009), charge
transfer (Côté, 2000), and cold chemistry, see (Pérez-Ríos, 2021) and
references therein.

Already at the level of the static properties, one encounters marked
differences due to the atom-ion interactions. The long-range attraction
between the neutral atoms and the ion leads to a number of new bound states,
which will be possible in the gas (Goold et al., 2010). This implies that
applicability of the standard mean-field theories such as the Gross-Pitaevski
approach must be re-evaluated (Massignan et al., 2005). At the level of
dynamics, the atom-ion systems present a number of opportunities as quantum
simulators of models of condensed-matter systems. As noted by (Bissbort et
al., 2013), a crystal of trapped ions immersed in a cold atomic gas features
close analogies to classic solid-state systems, and, thus, may provide a very
accurate quantum simulator for both structure and transport properties.
Furthermore, as recently pointed out, ions interacting with cold atoms in
tailored lattice setups may be used to simulate different types of lattice
gauge theories (Dehkharghani et al., 2017b) and Tomonaga-Luttinger liquids
(Michelsen et al., 2019).

In closing this section, we note that the combination of optical tweezers and
atom-ion systems has been recently proposed as a road to engineer spin-spin
interactions in more general quantum spin models (Arias Espinoza et al., 2021)
with possible applications in quantum simulation tasks or for quantum
computing purposes. For more information on different aspects of both
experimental and theoretical studies of atom-ion systems, we refer to the
review by (Tomza et al., 2019).

##  11 Summary and future prospectives

In this review we have discussed the physics of low-dimensional ultracold
quantum gases in the regime of few to many atoms, with the aim of providing a
comprehensive but still manageable account of the vast body of existing
literature. It is our hope that the review will be useful and serve as a guide
through the recent developments for both theoreticians and experimentalists
interested in this broad field of research and also inspire future endeavors.
Given that recent experimental advances represent a major motivation for the
research in this field, we have mainly focused on systems in harmonic traps
and pointed to the connections to homogeneous systems and lattice geometries
in the most relevant cases.

Many static properties of systems of interest are well understood by now, and
a large number of exact or approximate analytical methods exist that allow one
to build up some solid physical intuitions. In contrast, the corresponding
non-equilibrium quantum dynamics can often be analyzed only numerically, and
even in cases susceptible to exact treatments the inclusion of all
correlations is only possible within certain limits. This should, however, not
obscure the fact that recent progress in methodology has pushed the frontiers
enormously for systems in which all correlations can in principle be
calculated. Indeed, the work we have summarized in the review reveals that (i)
an in-depth analysis of the non-equilibrium dynamics in these controllable
quantum systems can deliver interesting insights and that (ii) ultracold atoms
may provide a highly valuable and flexible laboratory to study correlated
dynamical processes induced by a quench, continuous driving, or a suitable
time-dependent protocol. A comparison of experimentally obtained data with
state-of-the-art numerical results and advanced analytical models will keep
leading to a deeper understanding of non-equilibrium processes in the future,
including the physics of other interesting effects such as dissipation and
decoherence. It might even pave the way to prepare novel metastable states
with unseen properties or even functionalities.

Let us finally indicate a number of further research directions for many of
the systems discussed in the review. This is by no means an exhaustive list of
topics, but merely a selection mainly based on our own interests. It is our
humble attempt to create an outlook into the future from our limited
perspectives and the order in which we present our points should not be taken
as a prioritization.

To further understand the few- to many-body crossover, a new and promising
route is to investigate the impact of physical effects in spinor gases in
different correlation regimes. More generally, multicomponent systems in
different interaction regimes have recently been shown to hold great potential
to create up-to-now unseen quantum many-body states (Keller et al., 2022), for
which fundamental questions about the on-demand generation of spin-wave
excitations, triggering controllable spin-mixing processes and the impact of
correlations on the excitation spectrum exist. The latter has recently been
explored for many-body systems (Lavoine et al., 2021), and it would be
interesting to observe how these features will develop as one explores the
crossover down to few-body systems. Furthermore, in the context of spinor
gases, unraveling the emergence of static phases originating from the
cooperation of spin and spatial degrees-of-freedom as well as the competition
between intraband and interband tunneling dynamics in optical lattices are
topics worth pursuing.

Going beyond simple short-range interactions, we believe that it would be
fruitful to explore the impact of correlations on the interplay between
contact and long-range interactions during non-equilibrium dynamics. One
particular aspect here is how few-body mixtures can be used to design complex
pattern formation by triggering specific instabilities and excitation
mechanisms (Chomaz et al., 2022). More generally, the utilization of the
presence of quantum fluctuations in quantum liquid phases (droplets) as a
means to explore the fate of non-linear structures or the emergence of quasi-
particles ranging from polarons to magnons in systems where short- and long-
range interactions are competing is an intriguing problem.

Another challenging direction that has a long-standing history in both few-
and many-body physics, is the competition between quantum fluctuations and
temperature effects. From our point of view there is a need to further develop
novel effective models as well as numerical techniques in order to elucidate
this interplay. In conjunction with the additional degrees-of-freedom that a
multicomponent setting possesses, one may use few-body setups to elucidate new
applications within quantum information and quantum thermodynamics (Fogarty
and Busch, 2020), as well as the connection between information scrambling and
work statistics (Mikkelsen et al., 2022) to the spread of correlations in
small systems. Additionally, it would be worth exploring how issues of
topological physics in finite interacting systems can impact the physics and
the interplays discussed above, also in one-dimension (Reshodko et al., 2019).
Here we again stress that multicomponent setups are particularly interesting
as internal degrees-of-freedom can be reinterpreted as synthetic dimensions as
the particle number in the system grows. This poses the interesting question
of how this behavior evolves as a function of system size.

As we have discussed in previous Chapters, a key question for low-dimensional
systems is whether they are integrable or not. To further comprehend and
explore the transition from integrability to a non-integrable situation in
which chaos may arise, machine learning shows promise in identifying this
crossover (Kharkov et al., 2020; Huber et al., 2021). A major topic of
interest for this discussion is "how many interacting particles are needed to
bring a quantum system to the regime of strong chaos?" (Schiulaz et al., 2018;
Zisling et al., 2021; Fogarty et al., 2021). This question is an important
aspect to explore in few-body physics as it has implications on the study of
transport and quantum control, among others. It also has the potential to aid
the design and control quantum states and devices, and can open new paths as a
foundation for next-generation quantum technologies, through the use of
methods developed around the notion of quantum integrability and its breaking
(Wittmann Wilsmann et al., 2018; Grün et al., 2022).

However, above all, we must stress again the unstoppable and always prodigious
progress on the experimental side that will keep being a motivator for new and
advanced work also on the theory side. As one example, it is worth mentioning
the recently obtained control for allowing time-resolved detection of even
single atoms with high spatial resolution (Guan et al., 2019). These kind of
advanced techniques can open up new directions and shed new light also on old
results. They are part of the perpetual exchange of advances and ideas between
experimental and theoretical developments, that has been a constant and
fundamental pillar of supporting the fast progress in this area.

The advances mentioned above are some of the insights and results that we
believe may pave the way for future studies. These will address research
questions related to thermalization, quantum hydrodynamics, state engineering,
chaos and integrability and in general the management of correlations, and in
turn, these directions of research may become starting points for developments
in quantum technologies. We hope that our review will provide a useful
stepping stone for the corresponding theoretical developments and serve as a
fruitful guide for experimental advances.

## ACKNOWLEDGMENTS

This review could not have been written without the many fruitful discussions
and great collaborations with colleagues throughout the years, there are too
many to mention. Here we acknowledge conversations regarding the context of
the review with Lukas Rammelmüller, Oleksandr Marchukov, Panagiotis
Giannakeas, Tommaso Macrí and Manuel Valiente. S. I. M. acknowledges support
from the NSF through a grant for ITAMP at Harvard University. T.F.
acknowledges support from JSPS KAKENHI-21K13856 and T.F. and Th.B. acknowledge
support from the Okinawa Institute for Science and Technology Graduate
University. A.F. and R. E. B. acknowledge support from CNPq (Conselho Nacional
de Desenvolvimento Científico e Tecnológico). A. G. V. acknowledges support by
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie Grant Agreement No. 754411. P. S. is supported by the
Cluster of Excellence ’Advanced Imaging of Matter’ of the Deutsche
Forschungsgemeinschaft (DFG) - EXC2056 - project ID 390715994. N. T. Z. is
partially supported by the Independent Research Fund Denmark.

## References

  * Abad and Recati (2013) Abad, M., Recati, A., 2013.  A study of coherently coupled two-component Bose-Einstein condensates.  The European Physical Journal D 67, 1–11.  doi:10.1140/epjd/e2013-40053-2.
  * Abanin et al. (2019) Abanin, D.A., Altman, E., Bloch, I., Serbyn, M., 2019\.  Colloquium: Many-body localization, thermalization, and entanglement.  Rev. Mod. Phys. 91, 021001\.  doi:10.1103/RevModPhys.91.021001.
  * Agrawal (2000) Agrawal, G.P., 2000.  Nonlinear fiber optics, in: Nonlinear Science at the Dawn of the 21st Century. Springer, pp. 195–211.
  * Aikawa et al. (2012) Aikawa, K., Frisch, A., Mark, M., Baier, S., Rietzler, A., Grimm, R., Ferlaino, F., 2012.  Bose-Einstein Condensation of Erbium.  Phys. Rev. Lett. 108, 210401\.  doi:10.1103/PhysRevLett.108.210401.
  * Alba et al. (2019) Alba, V., Bertini, B., Fagotti, M., 2019.  Entanglement evolution and generalised hydrodynamics: interacting integrable systems.  SciPost Phys. 7, 5\.  doi:10.21468/SciPostPhys.7.1.005.
  * Albiez et al. (2005) Albiez, M., Gati, R., Fölling, J., Hunsmann, S., Cristiani, M., Oberthaler, M.K., 2005\.  Direct Observation of Tunneling and Nonlinear Self-Trapping in a Single Bosonic Josephson Junction.  Phys. Rev. Lett. 95, 010402\.  doi:10.1103/PhysRevLett.95.010402.
  * Alet and Laflorencie (2018) Alet, F., Laflorencie, N., 2018\.  Many-body localization: An introduction and selected topics.  Comptes Rendus Physique 19, 498–525.  doi:10.1016/j.crhy.2018.03.003.
  * Alon (2019) Alon, O.E., 2019.  Variance of a Trapped Bose-Einstein Condensate, in: J. Phys.: Conference Series, IOP Publishing. p. 012009.  doi:doi.org/10.1088/1742-6596/1206/1/012009.
  * Alon et al. (2006) Alon, O.E., Streltsov, A.I., Cederbaum, L.S., 2006.  Demixing of Bosonic Mixtures in Optical Lattices from Macroscopic to Microscopic Scales.  Phys. Rev. Lett. 97, 230403\.  doi:10.1103/PhysRevLett.97.230403.
  * Alon et al. (2008) Alon, O.E., Streltsov, A.I., Cederbaum, L.S., 2008.  Multiconfigurational time-dependent Hartree method for bosons: Many-body dynamics of bosonic systems.  Phys. Rev. A 77, 033613\.  doi:10.1103/PhysRevA.77.033613.
  * Alon, O. E. et al. (2004) Alon, O. E., Streltsov, A. I., Sakmann, K., Cederbaum, L. S., 2004\.  Continuous configuration-interaction for condensates in a ring.  Europhys. Lett. 67, 8–13.  doi:10.1209/epl/i2004-10047-3.
  * Altman et al. (2003) Altman, E., Hofstetter, W., Demler, E., Lukin, M.D., 2003\.  Phase diagram of two-component bosons on an optical lattice.  New J. Phys. 5, 113\.  doi:10.1088/1367-2630/5/1/113.
  * van Amerongen et al. (2008) van Amerongen, A.H., van Es, J.J.P., Wicke, P., Kheruntsyan, K.V., van Druten, N.J., 2008.  Yang-Yang Thermodynamics on an Atom Chip.  Phys. Rev. Lett. 100, 090402\.  doi:10.1103/PhysRevLett.100.090402.
  * Amico et al. (2017) Amico, L., Birkl, G., Boshier, M., Kwek, L.C., 2017\.  Focus on atomtronics-enabled quantum technologies.  New Journal of Physics 19, 020201\.  doi:10.1088/1367-2630/aa5a6d.
  * Amico et al. (2021) Amico, L., Boshier, M., Birkl, G., Minguzzi, A., Miniatura, C., Kwek, L.C., Aghamalyan, D., Ahufinger, V., Anderson, D., Andrei, N., Arnold, A.S., Baker, M., Bell, T.A., Bland, T., Brantut, J.P., Cassettari, D., Chetcuti, W.J., Chevy, F., Citro, R., De Palo, S., Dumke, R., Edwards, M., Folman, R., Fortagh, J., Gardiner, S.A., Garraway, B.M., Gauthier, G., Günther, A., Haug, T., Hufnagel, C., Keil, M., Ireland, P., Lebrat, M., Li, W., Longchambon, L., Mompart, J., Morsch, O., Naldesi, P., Neely, T.W., Olshanii, M., Orignac, E., Pandey, S., Pérez-Obiol, A., Perrin, H., Piroli, L., Polo, J., Pritchard, A.L., Proukakis, N.P., Rylands, C., Rubinsztein-Dunlop, H., Scazza, F., Stringari, S., Tosto, F., Trombettoni, A., Victorin, N., Klitzing, W.v., Wilkowski, D., Xhani, K., Yakimenko, A., 2021.  Roadmap on Atomtronics: State of the art and perspective.  doi:10.1116/5.0026178.
  * Anderegg et al. (2019) Anderegg, L., Cheuk, L.W., Bao, Y., Burchesky, S., Ketterle, W., Ni, K.K., Doyle, J.M., 2019.  An optical tweezer array of ultracold molecules.  Science 365, 1156–1158.  doi:10.1126/science.aax1265.
  * Andersen et al. (2016) Andersen, M.E.S., Dehkharghani, A.S., Volosniev, A.G., Lindgren, E.J., Zinner, N.T., 2016.  An interpolatory ansatz captures the physics of one-dimensional confined Fermi systems.  Scientific Reports 6, 28362\.  doi:10.1038/srep28362.
  * Andersen et al. (2017) Andersen, M.E.S., Harshman, N.L., Zinner, N.T., 2017.  Hybrid model of separable, zero-range, few-body interactions in one-dimensional harmonic traps.  Phys. Rev. A 96, 033616\.  doi:10.1103/PhysRevA.96.033616.
  * Andrade et al. (2021) Andrade, B., Kasper, V., Lewenstein, M., Weitenberg, C., Graß, T., 2021.  Preparation of the 1/2 Laughlin state with atoms in a rotating trap.  Phys. Rev. A 103, 063325\.  doi:10.1103/PhysRevA.103.063325.
  * Andriati and Gammal (2019) Andriati, A.V., Gammal, A., 2019\.  Superfluid fraction of few bosons in an annular geometry in the presence of a rotating weak link.  Phys. Rev. A 100, 063625\.  doi:10.1103/PhysRevA.100.063625.
  * Ao and Chui (1998) Ao, P., Chui, S.T., 1998.  Binary Bose-Einstein condensate mixtures in weakly and strongly segregated phases.  Phys. Rev. A 58, 4836–4840.  doi:10.1103/PhysRevA.58.4836.
  * Arias Espinoza et al. (2021) Arias Espinoza, J.D., Mazzanti, M., Fouka, K., Schüssler, R.X., Wu, Z., Corboz, P., Gerritsma, R., Safavi-Naini, A., 2021\.  Engineering spin-spin interactions with optical tweezers in trapped ions.  Phys. Rev. A 104, 013302\.  doi:10.1103/PhysRevA.104.013302.
  * Arkhipov et al. (2005) Arkhipov, A.S., Astrakharchik, G.E., Belikov, A.V., Lozovik, Y.E., 2005\.  Ground-state properties of a one-dimensional system of dipoles.  J. of Exp. and Th. Phys. Lett. 82, 39–43.
  * Armstrong et al. (2020) Armstrong, J.R., Jensen, A.S., Volosniev, A.G., Zinner, N.T., 2020\.  Clusters in separated tubes of tilted dipoles.  Mathematics 8.  doi:10.3390/math8040484.
  * Arwas et al. (2017) Arwas, G., Cohen, D., Hekking, F., Minguzzi, A., 2017\.  Resonant persistent currents for ultracold bosons on a lattice ring.  Phys. Rev. A 96, 063616\.  doi:10.1103/PhysRevA.96.063616.
  * Astrakharchik et al. (2008) Astrakharchik, G., Morigi, G., De Chiara, G., Boronat, J., 2008\.  Ground state of low-dimensional dipolar gases: Linear and zigzag chains.  Phys. Rev. A 78, 063622\.
  * Astrakharchik et al. (2021) Astrakharchik, G.E., Ardila, L.A.P., Schmidt, R., Jachymski, K., Negretti, A., 2021.  Ionic polaron in a Bose-Einstein condensate.  Communications Physics 4.  doi:10.1038/s42005-021-00597-1.
  * Astrakharchik et al. (2005) Astrakharchik, G.E., Boronat, J., Casulleras, J., Giorgini, S., 2005\.  Beyond the Tonks-Girardeau Gas: Strongly Correlated Regime in Quasi-One-Dimensional Bose Gases.  Phys. Rev. Lett. 95, 190407\.  doi:10.1103/PhysRevLett.95.190407.
  * Astrakharchik and Giorgini (2006) Astrakharchik, G.E., Giorgini, S., 2006\.  Correlation functions of a Lieb–Liniger Bose gas.  J. Phys. B: At. Mol. and Opt. Phys. 39, S1.  doi:10.1088/0953-4075/39/10/S01.
  * Astrakharchik and Malomed (2018) Astrakharchik, G.E., Malomed, B.A., 2018\.  Dynamics of one-dimensional quantum droplets.  Phys. Rev. A 98, 013631\.  doi:10.1103/PhysRevA.98.013631.
  * Atala et al. (2013) Atala, M., Aidelsburger, M., Barreiro, J.T., Abanin, D., Kitagawa, T., Demler, E., Bloch, I., 2013.  Direct measurement of the Zak phase in topological Bloch bands.  Nature Physics 9, 795–800.  doi:10.1038/nphys2790.
  * Atala et al. (2014) Atala, M., Aidelsburger, M., Lohse, M., Barreiro, J.T., Paredes, B., Bloch, I., 2014\.  Observation of chiral currents with ultracold atoms in bosonic ladders.  Nature Physics 10, 588–593.  doi:10.1038/nphys2998.
  * Baas et al. (2014) Baas, N.A., Fedorov, D.V., Jensen, A.S., Riisager, K., Volosniev, A.G., Zinner, N.T., 2014\.  Higher-Order Brunnian Structures and Possible Physical Realizations.  Physics of Atomic Nuclei. 77, 336–343.  doi:10.7868/S0044002714030027.
  * Baillie et al. (2017) Baillie, D., Wilson, R., Blakie, P., 2017.  Collective excitations of self-bound droplets of a dipolar quantum fluid.  Phys. Rev. Lett. 119, 255302\.
  * Bakr et al. (2009) Bakr, W.S., Gillen, J.I., Peng, A., Fölling, S., Greiner, M., 2009.  A quantum gas microscope for detecting single atoms in a Hubbard-regime optical lattice.  Nature 462, 74–77.
  * Banerjee (2007) Banerjee, A., 2007.  Collective oscillations of a Bose-Fermi mixture: Effect of unequal masses of Bose and Fermi particles.  Phys. Rev. A 76, 023611\.  doi:10.1103/PhysRevA.76.023611.
  * Baranov et al. (2012) Baranov, M., Dalmonte, M., Pupillo, G., Zoller, P., 2012\.  Condensed Matter Theory of Dipolar Quantum Gases.  Chem Rev. 112, 5012–5061.  doi:10.1021/cr2003568.
  * Barfknecht et al. (2015) Barfknecht, R.E., Brouzos, I., Foerster, A., 2015.  Contact and static structure factor for bosonic and fermionic mixtures.  Phys. Rev. A 91, 043640\.  doi:10.1103/PhysRevA.91.043640.
  * Barfknecht et al. (2016) Barfknecht, R.E., Dehkharghani, A.S., Foerster, A., Zinner, N.T., 2016\.  Correlation properties of a three-body bosonic mixture in a harmonic trap.  Journal of Physics B: Atomic Molecular Physics 49, 135301.  doi:10.1088/0953-4075/49/13/135301.
  * Barfknecht et al. (2017) Barfknecht, R.E., Foerster, A., Zinner, N.T., 2017.  Dynamical realization of magnetic states in a strongly interacting Bose mixture.  Phys. Rev. A 95, 023612\.  doi:10.1103/PhysRevA.95.023612.
  * Barfknecht et al. (2018a) Barfknecht, R.E., Foerster, A., Zinner, N.T., 2018a.  Effects of Interaction Imbalance in a Strongly Repulsive One-Dimensional Bose Gas.  Few-Body Systems 59, 22\.  doi:10.1007/s00601-018-1352-4.