the two-fold degenerate state $(0,4)$, as proposed in Ref. [73], is obtained
(see $\omega_{7}$ and $\omega_{8}$). The assignment of the corresponding
single GP state to $\omega_{7}$ is arbitrary, and one would need to analyse
for instance the response density to identify whether it correlates with the
7th or 8th many-body excitation. The overall accuracy of excitation energies
increases with the number of orbitals. Underlined digits denote deviations
from the exact values from Eqs. (160) and (161). All quantities are
dimensionless. See text for more details.

For the benchmark against the 1D repulsive HIM, the system parameters are set
to $\Omega=1$ for the trap frequency and $\lambda_{0}=-0.01$ for the
interaction strength, which yields a slight depletion of $f\approx 0.03\%$ for
$N=10$ bosons and $M=4$ orbitals. Calculations were carried out on a grid from
$[-9,9)$ with 128 grid points. Results for the excitation energies relative to
the ground-state energy, $\omega_{i}=E_{i}-\varepsilon^{0}$, are shown in
Table 2. There are several comments to be made. At first, it can be observed
that the mean-field BdG calculation ($M=1$) misses several excited states in
the low-energy spectrum. This is due to the fact that the BdG theory by
construction does not account for excitations were several particles at a time
are excited from the ground-state manifold, so-called multi-particle
excitations. For further details see Appendix A.1.1. The second observation is
that the numerical accuracy is enhanced if more self-consistent orbitals are
included into the ground-state description. Also the two-fold degenerate
excitation $(0,4)$, as proposed in Ref. [73], is obtained. The accuracy of the
BdG spectrum, apart from the c.m. excitation $(1,0)$, is clearly lower than
for the many-body results of $M=3$ and $4$ orbitals. In addition, the
convergence towards the exact ground-state energy is observed.

Table 3 shows excitation energies obtained for the isotropic 2D HIM with trap
frequencies $\Omega_{x}=\Omega_{y}=1$ and $N=100$ repulsive bosons with
interaction parameter $\lambda_{0}=-0.001$. All computations were carried out
on a $[-9,9)\times[-9,9)$ grid with $64\times 64$ grid points. At first, one
sees that the BdG approach yields the energy of the lowest c.m. excited state,
which is two-fold degenerate, to perfect accuracy. All other states obtained
are pure relative excitations. Higher c.m. excited states, as well as
combinations of those with relative excitations, are missing in the BdG
spectrum.

$\omega_{i}(\gamma)$ |  $M=1$ |  $M=2$ |  $M=3$ |  $(m_{x},m_{y},n_{x}$+$n_{y})$ | Exact
---|---|---|---|---|---
$\varepsilon^{0}(1)$ |  89.554453 |  89.551373 |  89.548293 | $(0,0,0)$ | 89.548292
$\omega_{1}(2)$ |  1.000000 |  1.000000 |  1.000000 | $(1,0,0)$ | 1.000000
|  1.000000 |  0.999792 |  1.000000 | $(0,1,0)$ | 1.000000
$\omega_{2}(3)$ |  1.791089 |  1.788859 |  1.788859 | $(0,0,2)$ | 1.788854
|  1.791089 |  1.791089 |  1.788861 | $(0,0,2)$ | 1.788854
|  1.791089 |  1.791143 |  1.788859 | $(0,0,2)$ | 1.788854
$\omega_{3}(3)$ |  n/a |  2.000438 |  2.000439 | $(2,0,0)$ | 2.000000
|  n/a |  n/a |  2.000657 | $(1,1,0)$ | 2.000000
|  n/a |  n/a |  2.000439 | $(0,2,0)$ | 2.000000
$\omega_{4}(4)$ |  2.686634 |  2.683292 |  2.683226 | $(0,0,3)$ | 2.683282
|  2.686634 |  2.684452 |  2.683226 | $(0,0,3)$ | 2.683282
|  2.686634 |  2.685519 |  2.683312 | $(0,0,3)$ | 2.683282
|  2.686634 |  2.686634 |  2.683312 | $(0,0,3)$ | 2.683282
$\omega_{5}(6)$ |  n/a |  2.789191 |  2.787088 | $(1,0,2)$ | 2.788854
|  n/a |  2.791314 |  2.787088 | $(1,0,2)$ | 2.788854
|  n/a |  2.793711 |  2.789270 | $(1,0,2)$ | 2.788854
|  n/a |  n/a |  2.789270 | $(0,1,2)$ | 2.788854
|  n/a |  n/a |  2.791536 | $(0,1,2)$ | 2.788854
|  n/a |  n/a |  2.791535 | $(0,1,2)$ | 2.788854
Table 3: Benchmark of the LR-MCTDHB implementation to the isotropic HIM in 2D.
Results are presented for $N=100$ bosons and different numbers of orbitals
$M$. The trapping frequencies are $\Omega_{x}=\Omega_{y}=1.0$, and the
interaction strength is $\lambda_{0}=-0.001$. The depletion is marginal with
$f\approx 0.003\%$ for $M=2$ and $f\approx 0.006\%$ for $M=3$. Shown are the
energies of the ground state $\varepsilon^{0}$ and the first $5$ excited
states, $\omega_{i}=E_{i}-\varepsilon^{0}$, together with their multiplicities
$\gamma$. In the BdG case, $M=1$, only the first pure c.m. excitation is
found, together with all pure excitations of relative coordinates. For LR-
MCTDHB(2), one obtains more c.m. excited states in the $x$-direction, together
with combinations of excitations in the relative coordinates. For LR-
MCTDHB(3), one obtains the missing c.m. excited states in the $y$-direction
and their combinations with excitations of the relative coordinates. The
overall accuracy increases with $M$, with few exceptions in the spectrum for
$M=2$ which are discussed in the text. Underlined digits denote deviations
from the exact values from Eqs. (160) and (161). All quantities are
dimensionless. See text for more details.

In contrast to the BdG approach, LR-MCTDHB(2) gives access to more excited
states. Before discussing those, it is worth examining the two underlying
ground-state orbitals. The first one is a Gaussian in the trap center, whereas
the second one resembles a $p_{x}$-orbital. The degree of depletion is
marginal with $f\approx 0.003\%$. In addition to the states obtained for
$M=1$, the spectrum contains the second-order c.m. excitation in the
$x$-direction, as well as three out of six states corresponding to
$\omega_{5}$, which are excited states combining a first-order c.m. excitation
in the $x$-direction and relative excitations. The accuracy of levels that
were already accessible within the BdG approach has increased, e.g., for the
first state of $\omega_{2}$ or the states of $\omega_{4}$. Nevertheless, some
states have lost accuracy. The most remarkable example is the second first-
order c.m. excitation of $\omega_{1}$, which is a c.m. excitation in the
$y$-direction. Naturally the question arises why this is the case. The reason
for this is that the second ground-state orbital, as discussed above,
resembles a $p_{x}$-orbital. This means that the ground-state description for
$M=2$ orbitals shows a preferred direction, which in this case is the
direction along the $x$-axis. Apparently, the consequence of this is that
excitations in the $y$-direction either lose in accuracy or remain
inaccessible. A possible reason for the loss in accuracy is that, compared to
the BdG case, a small amount of probability weight of the first orbital has
been transferred to the second orbital that effectively only covers the
$x$-direction. It is stressed that rotating the second orbital by an arbitrary
angle, e.g., $\pi/4$ or $\pi/2$, leads to exactly the same spectrum, which
means that it does not matter for the numerical results which direction is
preferred.

Including $M=3$ orbitals significantly enhances the obtained spectrum. One
observes that both problems of the spectrum for $M=2$, i.e., the loss in
accuracy for certain states as well as the absence of excitations in the low-
energy spectrum, are solved. The third ground-state orbital resembles a
$p_{y}$-orbital, i.e., it is perpendicular to the preferred direction of the
$M=2$ ground state. One can deduce from this observations a fundamental
difference between the computation of the ground state of a given $N$-boson
system and its excited states. Whereas increasing the number of orbitals
always leads to a lower ground state energy, it can happen that the spectrum
for ($M+1$) orbitals contains states that are less accurate than for only $M$
orbitals. This is in particular likely to happen in 2D and 3D because
preferred directions may occur. As a result, in order to obtain accurate
excitation spectra, one has to pay much more attention to the symmetry of the
problem in higher spatial dimensions.

To summarize, it was demonstrated for the HIM that both the MCTDHB and LR-
MCTDHB implementations utilized in this work are well-suited to obtain
accurate results for the spectrum and dynamics of trapped interacting bosons.
Furthermore, the many-body approaches clearly improve the GP and BdG mean-
field methods. As mentioned above, further benchmarks are given in Appendix C.

##  6 Applications to many-body dynamics and excitations of BECs

Since the advent of MCTDHB in 2007, the amount of scientific publications
where it is utilized grew substantially. As one of the first applications, the
splitting of 1D repulsive, fully-coherent BECs by a potential barrier has been
investigated [74, 75], showing that it can lead to fragmented condensates.
Similar observations were made by splitting a radially-symmetric BEC in 2D
with a ring-shaped barrier [76]. Also the reverse is possible, meaning that
two initially independent and thus fragmented BECs with no overlap in space
can collide, interfere and build up coherence [77]. Also for attractive BECs,
the splitting due to a Gaussian-shaped barrier results in the formation of a
superposition of two distinct bosonic clouds, a so-called caton, which is not
describable at the GP level [78]. Moreover, it has been shown for the first
time that coherent attractive BECs in 1D can fragment even without a potential
barrier, namely when the energy exceeds a threshold value [79]. Fragmentation
persists and may even intensify when more orbitals are included in the
calculation [80]. In the attractive case full convergence is difficult to
achieve as it requires more orbitals, see also [81]. When an attractive BEC
propagates towards a potential barrier and scatters from it, the system also
evolves to be highly fragmented although it was initially condensed [82, 83].
The degree of fragmentation depends on the degree of transmission through the
barrier. Even bright-soliton trains, i.e., stable multi-hump matter waves of
attractive bosons, whose dynamics were believed to be fully describable at the
mean-field level, are shown to lose their initial coherence [84].

The splitting of BECs has been further investigated with respect to the
generation of number- or phase-squeezed states applying optimal control theory
[85, 86, 87, 89, 88]. Optimal control has also been utilized in 1D bosonic
Josephson junctions (BJJ), on the one hand to steer the system parameters such
that an enhancement of the Shapiro effect was observed [90], and on the other
hand to drive a BEC from an initial state to a target state at the quantum
speed limit [91]. Moreover, it has been demonstrated that Mach-Zehnder
interferometry with ultracold bosons is relatively robust against the
nonlinear interaction between the particles [92].

A remarkable correspondence between the onset of wave chaos at the GP level
and the onset of fragmentation at the many-body level has been found for BECs
that scatter from either shallow periodic or disordered potential landscapes
[93, 94], but the result holds also for other potentials. Thus, the
development of wave chaos, i.e., the exponential separation of initially close
states, can be seen as an indication that a many-body treatment is necessary.
The distance between states is measured utilizing the $L^{2}$ norm between two
Hilbert space vectors.

Of potentially high relevance is the demonstration how the positions of
individual particles can be constructed from the many-body wave function by
simulating single shots [95]. The latter are connected to experiments on
trapped BECs where typically an image of the density represents a histogram of
a single shot.

Over the last years, MCTDHB has been applied to a growing number of systems,
e.g., to study the breathing dynamics in 1D harmonic traps due to a quantum
quench [96] or the many-body effects in solitonic excitations of ultracold
BECs [97, 98]. Other examples are the correlated dynamics of a single atom
coupled to a trapped BEC by collisions [99] or the examination of the
structure of mesoscopic molecular ions [100], which contributes to the
currently very active research in this field.

Recently, MCTDHB has also been applied to dipolar BECs in a double well [101]
and in 1D lattices [102, 103]. In addition, the coupling of one- and two-
component BECs in a cavity to a radiation field has been of interest [104,
105, 106]. One of the main results of the latter research was the finding of a
fragmented superradiant phase when the pump power of the laser exceeds a
critical threshold. This fragmented superradiance cannot be explained by the
Dicke model which assumes the BEC to be a simple two-level system.

In the remaining part of this section, the focus is laid on several
applications where the tunneling dynamics in traps (Section 6.1) and the
dynamical fragmentation of initially coherent BECs (Section 6.2) is
investigated, both in 1D and 2D systems. Afterwards, all applications of LR-
MCTDHB that exist up to now are presented in detail (Section 6.3). Special
emphasis for the 2D systems under consideration is laid on the impact of
angular momentum on the appearance of many-body effects in the dynamics and
excitation spectra.

###  6.1 Tunneling dynamics in traps

Subsequently, the tunneling dynamics of BECs held in double-well traps in both
one and two spatial dimensions are discussed in this section. The main purpose
is to apply MCTDHB to solve the many-boson Schrödinger equation and to
discover novel phenomena that are not captured by the standard methods.

####  6.1.1 Exact tunneling dynamics in a bosonic Josephson junction

The quantum dynamics of a tunneling BEC in a 1D BJJ have constituted a large
research field in the past decades. Before the first experimental realizations
[107, 108], theoretical predictions concerning Josephson oscillations and
self-trapping, i.e., the suppression of tunneling between the wells, were
obtained by a (multi-mode) GP description [109, 110, 111, 112, 113, 114], a
path-integral formulation [115], a two-mode Bose-Hubbard (BH) approach [116]
or a Fock-space WKB method [117]. A review of various theoretical and
experimental results is given in [118]. A more recent study dealt with the
tunneling dynamics in a controlled BJJ, where the spin state of an ionic
impurity controls the tunneling between the wells [119]. The latter extends
the findings on the dynamics of trapped ultracold BECs due to the presence of
an ion [120, 121]. Furthermore, the orbital Josephson effect due to a time-
dependent driving potential was investigated numerically by employing MCTDHB,
not only in a double-well system [122].

The main objective of this section is however to study the out-of-equilibrium
tunneling dynamics of a single-component BEC in a BJJ from a many-body
perspective, and the method of choice is MCTDHB. As described below, the
latter approach allows for the observation of fundamentally new physics
dealing with, e.g., reduced self-trapping and dynamical fragmentation of
initially fully-coherent bosonic clouds. A comparison to other popular
theoretical tools like the GP approach and the BH model is made. The findings
presented below were recently published in Refs. [123, 124, 125, 126], where
further details can be found. Recently, results of a comparable study where
both the BH model and MCTDHB were used to analyze the tunneling of few-boson
systems in asymmetric double-wells, investigating the interaction blockade
that isolates the motion of a single particle in the vicinity of others,
became available [127].

The trapping confinement considered here reflects two harmonic wells of the
form

$V_{\pm}(x)=\frac{1}{2}(x\pm 2)^{2}$ (172)

which are connected by a cubic spline in the interval $|x|\leq 0.5$ that
generates a potential barrier between the wells. The plus sign refers to the
left well, whereas the minus sign refers to the right well. The repulsion
between the bosons is modeled by the contact interaction potential,
$\lambda_{0}\hat{W}(x,x^{\prime})=\lambda_{0}\delta(x-x^{\prime})$. The
repulsion strength is given by the mean-field parameter
$\Lambda=\lambda_{0}(N-1)$. An important quantity with respect to the
tunneling dynamics is the occupation probability of the left and right wells,
where the former is given by

$P_{L}(t)=\frac{1}{N}\int_{-\infty}^{0}\rho(x,t)dx$ (173)

where $\rho(x,t)$ describes the time-dependent one-body density, see Eq.
(2.4), and $N$ refers to the number of particles in the trap.

Since the MCTDHB results are compared to the BH model below, it is instructive
to define the left- and right-well orbitals, $\phi_{L}$ and $\phi_{R}$, with
the ground and the first excited state of the trap, denoted by $\phi_{g}$ and
$\phi_{u}$ due to their gerade and ungerade symmetry. These two states are
energetically well below the barrier. The orbitals $\phi_{L}$ and $\phi_{R}$
are defined via

$\phi_{L,R}(x)=\frac{\phi_{g}(x)\pm\phi_{u}(x)}{\sqrt{2}}.$ (174)

The BH parameters for on-site interaction $U$ and hopping $J$ are derived from
$\phi_{L}$ and $\phi_{R}$ by

$U=\lambda_{0}\int|\phi_{L}(x)|^{4}\,dx,\quad
J=-\int\phi_{L}^{\ast}(x)\hat{h}(x)\phi_{R}(x)\,dx$ (175)

and utilized to define the BH interaction parameter
$\Lambda_{\text{BH}}=UN/(2J)$. The two-site BH Hamiltonian then reads

$\hat{H}_{\text{BH}}=-J(\hat{b}_{L}^{\dagger}\hat{b}_{R}+\hat{b}_{R}^{\dagger}\hat{b}_{L})+\frac{U}{2}(\hat{b}_{L}^{\dagger}\hat{b}_{L}^{\dagger}\hat{b}_{L}\hat{b}_{L}+\hat{b}_{R}^{\dagger}\hat{b}_{R}^{\dagger}\hat{b}_{R}\hat{b}_{R})$
(176)

where the operators $\hat{b}_{L,R}^{(\dagger)}$ annihilate (create) a boson in
the left and right well, respectively. Within the two-mode BH approach, the
eigenvalues and eigenfunctions of the one-body RDM

$\rho^{(1)}=\begin{pmatrix}\rho_{LL}&\rho_{LR}\\\
\rho_{RL}&\rho_{RR}\end{pmatrix}$ (177)

with $\rho_{LL}=\langle
N,0|e^{+i\hat{H}_{\text{BH}}t}\,\hat{b}^{\dagger}_{L}\hat{b}_{L}\,e^{-i\hat{H}_{\text{BH}}t}|N,0\rangle$
(with $\hbar=1$) and the other elements defined analogously, represent the
natural occupations and natural orbitals of the BH method. Therefore, the
population of the second natural orbital can be used as a measure for
fragmentation.

In the following, the tunneling dynamics when the BEC is initially trapped in
the left well, i.e., $P_{L}(0)=1$, will be explored for different repulsion
strengths. Most importantly, for a two-mode GP description of the system,
self-trapping of the bosons is predicted for
$\Lambda_{\text{BH}}>\Lambda_{c}=2$ when the magnitude of the initial
imbalance between the occupations of the two wells, denoted by
$z=\frac{N_{L}-N_{R}}{N}$, is unity [109, 110, 111].

Figure 5: Short-time dynamics of the occupation probability of the left well,
$P_{L}(t)$, for (a) $N=20$, $\Lambda=0.152$, (b) $N=100$, $\Lambda=0.152$
($\Lambda_{\text{BH}}<\Lambda_{c}$) (c) $N=20$, $\Lambda=0.245$, and (d)
$N=100$, $\Lambda=0.245$ ($\Lambda_{\text{BH}}>\Lambda_{c}$). The GP [dotted
black], BH [dashed magenta] and numerically exact results from MCTDHB($M$)
[solid blue] are compared. For the weak interaction in panels (a) and (b), the
GP theory predicts Rabi oscillations of the full density. At the many-body
level, a collapse of these oscillations occurs after approximately three
cycles. For stronger interaction in panels (c) and (d), GP still predicts that
all bosons tunnel between the wells, whereas at the many-body level tunneling
is clearly suppressed. However, MCTDHB and the BH model yield different
dynamics in all cases. The deviations become more pronounced for stronger
repulsion. The inset shows the numerical convergence with respect to $M$ for
2, 4, 6 and 8 orbitals [(a) and (c)] and 2 and 4 orbitals [(b) and (d)]. All
quantities are dimensionless. See text for details. The figure is taken from
Ref. [123]. Figure 6: Upper panel: Dynamics of $P_{L}(t)$ for very strong
repulsion. The interaction parameter $\Lambda_{\text{BH}}=47.8$ ($43.4$) for
$N=10$ [solid blue] ($N=100$, solid green) is clearly above the critical value
of $\Lambda_{c}=2$. While the BH predicts full self-trapping of the bosons
[dashed magenta], the exact results show tunneling between the wells. Lower
panel: Evolution of the natural occupations $n_{i}^{(1)}(t)$. At the BH level,
the system remains condensed throughout the propagation time, while it
fragments at the exact many-body level where essentially four orbitals become
macroscopically occupied. All quantities are dimensionless. See text for
details. The figure is adapted from Ref. [123]. Figure 7: BH versus exact
dynamics for both attractive and repulsive BECs with system parameters $N=20$,
$\frac{|U|}{J}=0.226$ and $|\lambda_{0}|=0.0129$. Left panel: Occupation
probability $P_{L}(t)$ of the left well. While the BH predicts the same
evolution for attraction and repulsion, differences can be observed at the
exact many-body level. Inset: Time-evolution of the ratio of $P_{L}(t)$ for
attraction and repulsion. Right panel: The same for the first two natural
occupation numbers $n_{i}^{(1)}(t)$. Again, the BH model predicts equivalent
dynamics for attraction and repulsion, whereas the exact results show clear
differences. All quantities are dimensionless. See text for details. The
figures are taken from Ref. [124].

Fig. 5 shows a comparison between the short-time dynamics of $P_{L}(t)$ for
BECs with $N=20$ ($100$) bosons obtained from MCTDHB (solid blue), BH (dashed
magenta) and GP (dotted black) for two different interaction strengths where
one is below and one is above the critical value $\Lambda_{c}$. One observes
that for the case of $20$ bosons with weak interaction parameter
$\Lambda_{\text{BH}}<\Lambda_{c}$ [panel (a)], the GP prediction is totally
different than the BH and numerically exact MCTDHB predictions. According to
GP, Rabi oscillations of the entire cloud between the wells occur, whereas the
BH model and MCTDHB predict a density collapse after approximately three Rabi
cycles, resulting in roughly $50\%$ of the bosons in each well. The time for
one Rabi cycle is close to the theoretical prediction of
$t_{\text{Rabi}}=\pi/J$. Although the BH curve is clearly superior to the GP
curve, it also deviates from the numerically exact result after half a Rabi
cycle. In the case of $100$ bosons [panel (b)], similar observations can be
made. However, the density collapse occurs after a longer time.

As mentioned above, for stronger repulsion with
$\Lambda_{\text{BH}}>\Lambda_{c}$, tunneling should be entirely suppressed due
to self-trapping according to a two-mode GP description. It can be deduced
from the exact dynamics of $20$ bosons [panel (c)] that indeed less bosons
tunnel between the wells, meaning that the oscillation amplitude is clearly
smaller than for the weaker repulsion. As predicted by MCTDHB, at most $50$
bosons tunnel into the right well after half a Rabi cycle, and even less in
the BH dynamics. Nevertheless, tunneling is not completely suppressed. This
means that the effect of self-trapping, as predicted by the two-mode GP
theory, is obviously reduced at the accurate many-body level. The deviations
between the MCTDHB and BH results are more pronounced for stronger than for
weaker repulsion. The dynamics of $N=100$ bosons yields similar observations
[panel (d)]. In contrast to that, the GP theory predicts tunneling of the
whole cloud.

The reason why GP fails to describe the dynamics correctly is the development
of fragmentation, already on the short time scales shown. For the case of weak
repulsion in the upper panels of Fig. 5, the initially condensed BEC with
$f=10^{-4}$ ($10^{-5}$) for $N=20$ ($100$) evolves to be fragmented by $33\%$
($26\%$) after three Rabi cycles. Interestingly, for stronger repulsion, the
degree of fragmentation of the again initially fully-coherent BEC is less
after three Rabi cycles than for the case of weaker repulsion, yielding
approximately $28\%$ and $18\%$ for $20$ and $100$ bosons, respectively. The
predicted GP dynamics differ both quantitatively and qualitatively from the
exact results. Most important, for all system parameters used in Fig. 5, GP
predicts that essentially all bosons tunnel back and forth between the wells.
This is clearly not the case at the many-body level, neither for BH nor for
MCTDHB. It is stressed that although the chosen system parameters are within
the expected regime of validity for both the GP and BH approaches, both
theories fail to describe the short-time dynamics of the BECs.

To further illustrate the effect of reduced self-trapping as well as the
degree of fragmentation, Fig. 6 shows the evolution of $P_{L}(t)$ and of the
occupation of the first few natural orbitals $n_{i}^{(1)}(t)$, both for MCTDHB
and the BH model. The utilized repulsion parameter exceeds the critical value
significantly with $\Lambda_{\text{BH}}=47.8\,(43.4)$ for $N=10\,(100)$. From
the top panel, it can be inferred that the BH model predicts complete self-
trapping for both values of $N$, meaning that all bosons remain in the left
well throughout the time-evolution. On the contrary, the predictions of MCTDHB
show that indeed bosons tunnel from left to right, and one can anticipate that
$P_{L}(t)$ approaches the long-time average of $0.5$. With respect to the
natural occupations, one observes that the BEC remains fully-condensed in time
at the BH level, while the MCTDHB results show that fragmentation sets in
already during the first Rabi cycle. Moreover, more than ten orbitals are
necessary to accurately describe the dynamics. The different results are
associated with a quick loss of coherence between the bosons that is only
described by MCTDHB and not at the BH level. In this context, see Ref. [123]
for more details.

Another interesting aspect is the symmetry of the two-mode BH Hamiltonian in
Eq. (176) under the unitary transformation
$\hat{R}=\\{\hat{b}_{L}\rightarrow\hat{b}_{L},\hat{b}_{R}\rightarrow-\hat{b}_{R}\\}$,
leading to $\hat{R}\hat{H}_{\text{BH}}(U)\hat{R}=-\hat{H}_{\text{BH}}(-U)$. A
consequence of this is the equivalence of the occupation of the wells for
attractive and repulsive interactions, i.e.,

$P_{L,R}(t;U)=P_{L,R}(t;-U).$ (178)

Moreover, it can be shown that the eigenvalues of the BH one-body RDM in Eq.
(177) do not depend on the sign of the interaction. Thus, $P_{L}(t)$ and the
occupations $n_{i}^{(1)}(t)$ will be the same for attraction and repulsion of
the same magnitude. Naturally, the full many-body Hamiltonian of the system
does not possess such a symmetry, and one can therefore expect different
dynamics for attraction and repulsion between the bosons.

The left panel of Fig. 7 shows the short-time dynamics of $p_{L}(t)$ for a BEC
with $N=20$ bosons with interaction parameter $|U|/J=0.226$, corresponding to
$|\lambda_{0}|=0.0129$. As anticipated, the evolution at the BH level is
identical for attraction and repulsion. However, the results obtained by
MCTDHB show differences between the cases of positive and negative
$\lambda_{0}$. On the time scale shown, the BH model mainly underestimates
tunneling for the repulsive case, whereas it overestimates tunneling for the
attractive case. The right panel of Fig. 7 shows the evolution of the first
few natural occupation numbers for the same system as in the left panel, both
at the BH and full many-body level. As expected, the BH model predicts the
same evolution for $n_{1}^{(1)}(t)$ and $n_{2}^{(1)}(t)$, irrespective of the
sign of the interaction parameter. In contrast to that, the results obtained
from MCTDHB show that the evolution of the different occupations are distinct
for the attractive and repulsive cases. Compared to the exact results, the BH
model mainly overestimates the fragmentation of the attractive BEC, whereas it
mainly underestimates it for the repulsive BEC.

Figure 8: Universality of the degree of fragmentation. Upper panel: Evolution
of $P_{L}(t)$ for different numbers of bosons ($N=100$ in blue, $N=1000$ in
red and $N=10000$ in orange) where the interaction parameter
$\Lambda=\lambda_{0}(N-1)=0.152$ is kept fixed. Intermediate particle numbers
($N=200,\,500,\,2000$ and $5000$) are shown in gray. The less particles are
contained in the BEC, the faster the density collapses. Lower panel: Evolution
of the first two natural occupations. Irrespective of the number of bosons,
the degree of fragmentation upon the density collapse is the same. All
quantities are dimensionless. See text for details. The figure is taken from
[126].

The analysis of the dynamics in a 1D BJJ is closed by discussing the
compelling feature of universal fragmentation. The latter means that all BECs
with identical interaction parameter $\Lambda$ fragment to the same value
after the density has collapsed. Fig. 8 shows exactly this behavior for
$\Lambda=0.152$ and BECs with $N=100$, $1000$, and $10000$ particles where
initially all bosons are kept in the left well, i.e., $P_{L}(0)=1$. From the
upper panel, one sees that the density oscillations between the wells collapse
after a few Rabi cycles, leading to the equal distribution of $50\%$ of the
bosons in each well. In general, the less bosons are contained in the BEC, the
faster the density collapses. Coming along with the damping of the
oscillations, the BEC starts to fragment, reaching its maximal value after the
collapse. Again, the less bosons are in the system, the faster the system
fragments. Most importantly, irrespective of $N$, the final degree of
fragmentation is the same. The fact that the degree of fragmentation appears
to be universal is unexpected, since usually fragmentation strongly depends on
the number of particles in the BEC. Moreover, the degree of fragmentation upon
the density collapse, denoted by $f_{\text{col}}$, depends highly on the
initial imbalance $z$ of bosons in the left and right wells. The larger $z$,
the stronger does the system fragment in time. This holds true even in the
limit of very weak interactions, $\Lambda\ll 1$, meaning that there is no
weak-interaction regime where GP yields the correct dynamics. It is important
to note that the universality is not an artifact of any approximation, since
the result is obtained by solving the full many-boson Schrödinger equation
numerically exact. Interestingly, the two-mode BH approach can be used to
derive an analytic prediction of $f_{\text{col}}$ yielding results close to
the exact ones obtained by MCTDHB. Additional details on the analytic
expression can be found in Ref. [126].

To summarize, it was demonstrated that the out-of-equilibrium tunneling
dynamics of BECs in a 1D BJJ show significant many-body features that are not
accounted for at the GP mean-field level. This includes the appearance of a
density collapse in the oscillations between the wells as well as dynamical
fragmentation which the GP theory cannot capture by construction. The degree
of fragmentation is universal for fixed $\Lambda=\lambda_{0}(N-1)$ for a wide
range of initial states. Applying the commonly used BH model, it was shown
that it is superior in certain aspects in comparison to the GP mean-field
approach. However, it does also not predict the numerically exact tunneling
dynamics as obtained by MCTDHB, e.g., on the short-time dynamics where the
interaction parameter is in the vicinity of its critical value or for very
strong repulsion where it predicts complete self-trapping of the bosonic
cloud. However, the BH model can be utilized to derive an analytic expression
for the universal degree of fragmentation that is in good agreement with the
exact results.

####  6.1.2 Many-body tunneling dynamics in a two-dimensional radial double
well

After having shown the many-body nature of the tunneling dynamics of BECs in a
1D BJJ, the question of how a repulsive BEC behaves in a radially symmetric 2D
system with a double-well structure is addressed. In particular, as for the
previous section, the out-of-equilibrium tunneling dynamics are of interest,
and it is studied whether many-body effects can be observed. The results
presented in this section were published in Ref. [128], and additional details
can be found therein.

Figure 9: (a) Schematic plot of the radial double well. The IN region denotes
the trap center, the OUT region denotes the outer annulus. Both are separated
by the central ring-shaped barrier. (b) Ground-state energy per particle of a
BEC with $N=100$ bosons and total angular momentum $L=0$ for different
positions $R_{B}$ of the radial barrier. The results are calculated at the
mean-field GP level ($M=1$). The quantities $E_{0}^{\text{IN}}$ (short-dashed
red) and $E_{0}^{\text{OUT}}$ (solid green) denote the energies for the IN and
OUT regions for interaction strength $\Lambda=0$, whereas $E_{2}^{\text{IN}}$
(dash-double-dotted blue) and $E_{2}^{\text{OUT}}$ (dash-dotted magenta)
denote the energies for $\Lambda=2$. The crossing points, i.e., the radii at
which the energies for IN and OUT are alike, are located at $R_{1}=3.271$
($\Lambda=0$) and $R_{2}=3.412$ ($\Lambda=2$), see the vertical black lines.
Larger repulsion between the bosons shifts the crossing point to larger radii.
(c) Same analysis as in (b) but for total angular momentum $L=N$, meaning that
each particle carries on average one unit of angular momentum. The location of
the crossing points are $R_{3}=4.089$ ($\Lambda=0$) and $R_{4}=4.093$
($\Lambda=0.2$). The angular momentum has a significant effect on the position
of the crossing points. Increasing $L$ shifts them to larger radii, which is
qualitatively similar as for increasing $\Lambda$. However, one would need to
highly increase the interaction strength in order to shift the crossing points
as much as adding one unit of angular momentum per particle does. All
quantities are dimensionless. See text for details. Figs. (b) and (c) are
adapted from Ref. [128].

To study tunneling phenomena in 2D geometries has been of high interest in
recent years, especially for trapped vortices, i.e., macroscopic quantum
states with definite angular momentum, characterized by a density node at its
core and phase discontinuities around it. Research concerning tunneling of
vortices has been carried out for example in 2D superfluids [129]. Others
dealt with vortices tunneling through a Gaussian-shaped barrier [130], in
pinning potentials [131, 132, 133], or between two Gaussian wells [134], to
mention only a few. However, there has been no investigation of the tunneling
dynamics of vortices in a radially-symmetric double well, in particular not at
the many-body level. The latter geometry is very interesting to study many-
body dynamics since the trap symmetry preserves the angular momentum of the
bosonic cloud. It has been utilized recently to study many-body vortices which
are separated in space due to the central barrier [135]. The following
results, recently published in Ref. [128], especially deal with the impact of
angular momentum on the tunneling dynamics.

Accurate many-body dynamics are obtained by solving the full many-body
Schrödinger equation utilizing MCTDHB and compared to the corresponding mean-
field dynamics from the GP theory.

For the dynamics, the radial double well in which a BEC with $N=100$ spinless
bosons is confined reads

$V(r)=\begin{cases}B\,e^{-2(r-R_{B})^{4}}+C\,e^{-0.5\,(r-R_{C})^{4}}\,\,&\mbox{if
}r=(x^{2}+y^{2})^{1/2}\leq R_{C}\\\ C\,\,&\mbox{if }r>R_{C}\end{cases}$ (179)

where $B$ is the height of the central ring-shaped barrier located at the
radius $R_{B}$, and $C$ is the height of the crater wall at radius $R_{C}$. A
schematic illustration of $V$ is shown in Fig. 9(a). The region enclosed by
the barrier is denoted by IN, the external rim, i.e., the region between the
barrier and the crater wall, is denoted by OUT. For the numerical
calculations, $C=200.0$ and $B=1.0$ are set which ensures on the one hand that
the energy per particle is sufficiently small such that the BEC stays trapped
in the crater, and on the other hand allows for tunneling between IN and OUT,
as described below. The radius of the crater is set to $R_{C}=9.0$.

The short-range repulsion between the particles is modeled by a Gaussian that
has been employed recently [76, 136],

$\lambda_{0}\hat{W}(\mathbf{r}_{i},\mathbf{r}_{j})=\frac{\lambda_{0}}{2\pi\sigma^{2}}\,e^{-|\mathbf{r}_{i}-\mathbf{r}_{j}|^{2}/2\sigma^{2}},$
(180)

which circumvents the regularization problems of the contact interaction
[137]. The repulsion strength will again be expressed in terms of the mean-
field parameter $\Lambda=\lambda_{0}(N-1)$ in the following. The effect of the
range of the interaction on the degree of fragmentation has been investigated
in detail for repulsive BECs in 1D and 2D single wells [138], essentially
showing that it becomes increasingly dependent on the density the more long-
ranged it is.

The analysis starts with the ground-state energy of two distinct cases, namely
for (i) having all particles in the IN region with closed external rim and
(ii) having all particles in the OUT region with closed trap center. The
notation $E^{\text{IN}}_{\Lambda}$ and $E^{\text{OUT}}_{\Lambda}$, where the
superscript denotes whether the BEC is held in the IN or OUT region and the
subscript denotes the repulsion strength $\Lambda$, is utilized. The ground-
state energies for different radial positions of the barrier $R_{B}$ are
calculated using the GP equation. The corresponding results are shown in Fig.
9(b) for total angular momentum $L=0$ and Fig. 9(c) for $L=N$, i.e., where
each boson on average carries one unit of angular momentum. One observes that
for all shown parameters ($\Lambda=0,\,2.0$ for $L=0$ and $\Lambda=0,\,0.2$
for $L=N$) the curves for $E^{\text{IN}}_{\Lambda}$ and
$E^{\text{OUT}}_{\Lambda}$ intersect at certain radii which are referred to as
the crossing points or simply the crossings in the following. These are
defined as the positions where the ground-state energies of the BEC in the IN
and OUT regions are identical. Putting the barrier onto the crossing point is
the 2D analogue of the symmetric double well in 1D. However, the IN and OUT
regions are only energetically equivalent, whereas their geometries are
totally different. It is found that both the interaction strength $\Lambda$
and the angular momentum $L$ push the crossing point to larger radii once they
are increased. Though, the effect of $L$ is significantly stronger since
already for $\Lambda=0$, the transition from $L=0$ to $L/N=1$ leads to a shift
of the crossing point from $R_{B}=3.271$ to $R_{B}=4.089$. In order to obtain
a similar shift just by increasing $\Lambda$, the repulsion needs to be
largely increased, and thus one would certainly leave the weakly-interacting
regime. The depletion of the ground states with $M=4$ orbitals are $f\approx
10^{-3}$ for $L=0$ and $\Lambda=2.0$ and $f\approx 10^{-9}$ for $L=N$ and
$\Lambda=0.2$, respectively. This means that in both cases, the systems’
ground states are essentially condensed.

(a)

(b)

(c)

(d)

Figure 10: (a) Time evolution of the occupation probability of the external
rim, $P_{\text{OUT}}(t)$, for $N=100$ bosons with $L=0$ and $\Lambda=2.0$. For
$M=4$, initially almost all bosons tunnel. The oscillation amplitude is damped
over time and effectively collapses after 12 Rabi cycles with approximately
half of the bosons in each part (solid red). The system evolves from being
almost fully condensed to two-fold fragmented (see solid green and dark blue
curves for $n_{1}$ and $n_{2}$). The occupations of the third and fourth
orbital are identical and negligible (solid light blue). On the contrary, the
GP theory predicts undamped oscillations (dotted gray). (b) Same as in (a) but
for $L=N$ and $\Lambda=0.2$. Again, the system evolves from being fully
coherent to two-fold fragmented, and the particles finally occupy both the IN
and OUT parts by roughly 50$\%$. Mind the different time scales
($\tau_{1}>\tau_{2}$). Inset: Detailed view between $150\tau_{2}\leq t\leq
155\tau_{2}$. The GP dynamics again do not show the density collapse. (c)
Real, imaginary, and absolute value of the first two natural orbitals
$\alpha_{1}$ and $\alpha_{2}$ at $t=18.4\,\tau_{1}$ where the density is
already collapsed. The GP case is shown in the top panels. For $\alpha_{1}$
and $\alpha_{2}$, the density occupies both the IN and OUT regions. At the
same time, the GP orbital only covers the external part. (d) Same as in (c)
but for the dynamics described in (b). At the many-body level, the nodal
structure of the vortex state is clearly visible in the trap center. All
quantities are dimensionless. See text for details. All figures are adapted
from Ref. [128].

To study the out-of-equilibrium dynamics for different values of $\Lambda$ and
$L$, the focus is laid on the situation where the barrier is located on the
corresponding crossing points. First, the case of zero angular momentum, i.e.,
$L=0$, is analyzed. As an initial state, the ground state of trapping the BEC
in the OUT region with closed central part is chosen. To trigger the dynamics,
the system is quenched by suddenly opening the trap center, such that the
entire potential is given by Eq. (179). It was found numerically that there is
no difference in starting from either the IN or OUT region when the barrier is
located at the crossing point. Fig. 10(a) shows the dynamical evolution of the
occupation probability of the external rim,

$P_{\text{OUT}}(t)=1-P_{\text{IN}}(t)=\frac{1}{N}\int_{R_{B}<r\leq
R_{C}}\rho(\mathbf{r};t)d\mathbf{r},$ (181)

where $\rho(\mathbf{r},t)$ denotes the single-particle density. For the first
two Rabi cycles of $\tau_{1}=110.23$, nearly the whole cloud tunnels between
IN and OUT. Then, the occupation of the second natural orbital starts to grow,
accompanied by a damping of the tunneling oscillations. After approximately 12
cycles, the density oscillations have collapsed and the particles are almost
equally distributed between IN and OUT. The BEC is now highly (two-fold)
fragmented, with occupations of 56.7% and 41.1% of the first two natural
orbitals $\alpha_{1}$ and $\alpha_{2}$, respectively. In contrast to that, the
GP dynamics do not show any damping of the oscillations, not even for the
long-time behavior (not shown). Fig. 10(c) shows the real and imaginary parts
as well as the density of the GP orbital (top panels) and of the first two
natural orbitals in the many-body case (middle and lower panels) at an instant
in time at which the density oscillations have already collapsed. Whereas
there are no bosons in the trap center for the GP case, there is significant
population of the IN region at the many-body level.

Similar observations can be made for $L=N$ and $\Lambda=0.2$, presented in
Fig. 10(b). The many-body dynamics show that the density collapses and that
approximately 50 bosons are occupying the IN and OUT regions finally. The BEC
is two-fold fragmented, the occupation of higher natural orbitals is
negligible. The GP dynamics (see inset) do not account for the damped
oscillations. With regard to the density of the GP and natural orbitals in
Fig. 10(d), it can be seen that $\alpha_{1}$ and $\alpha_{2}$ clearly show the
characteristic nodal vortex structure in the center, whereas at the GP level
there are no bosons in the IN region at the chosen instant of time
$t=183.3\,\tau_{2}$ with $\tau_{2}=66.79$.

It was further found that many-body effects like the density collapse or the
onset of fragmentation are more likely to occur when the particles carry
angular momentum. For $L=N$ and $\Lambda=2.0$, the system quickly becomes
four-fold fragmented, i.e., at least four orbitals are necessary to account
for the many-body dynamics. For $L=0$, the dynamics are essentially governed
by only two orbitals, as can be seen from Fig. 10(c). Thus, the weakly-
interacting regime is apparently smaller for the case of $L>0$, meaning that
vortices are more sensitive to many-body effects. In Section 6.3.3, this will
be discussed again when excitations in rotating BECs are investigated.

To summarize, the many-body nature of the tunneling dynamics of BECs in a 2D
radial double well was demonstrated . The GP mean-field theory does not
account for the collapse of the density oscillations between IN and OUT, and
therefore fails to give an accurate description of the dynamics. Furthermore,
the initially condensed BECs fragment. It was found that bosons carrying
angular momentum are more sensitive to many-body effects than non-rotating
ones because fragmentation sets in already for weaker interaction strengths.
In agreement with these findings, a recent work demonstrates analytically that
a many-body approach is necessary for describing the conservation of angular
momentum in a 2D radially-symmetric system [139].

###  6.2 Dynamical fragmentation

####  6.2.1 Tunneling to open space

The obtained findings on the many-body tunneling process of initially trapped
bosons to open space in 1D are summarized in this section. The tunneling
process of a single particle was understood already at the beginning of the
20th century [140, 141, 142, 143] and is still of current interest, e.g., with
respect to the exit time and momentum of an electron in the ionization process
[144, 145]. Tunneling of a many-boson system, which has also been studied
[146, 147, 148, 149], is naturally more complicated due to the interactions
and correlations between the constituent particles. The following results,
published in Refs. [150, 151, 152, 153, 154], explain the tunneling mechanism
to open space from a many-body point-of-view utilizing MCTDHB. The main points
that are addressed deal with the questions whether the particles tunnel one by
one or if several particles tunnel simultaneously, and whether the entire
process is indeed of many-body nature.

To study the dynamics of the tunneling to open space, a system of $N$ bosons
confined in a harmonic trap of the form $V(x,t=0)=\frac{1}{2}x^{2}$ is
considered. For $t>0$, the trap is quenched such that it becomes open to one
side. The region where the bosons are trapped initially and the open space are
separated by a potential barrier. The full potential reads

$V(x,t>0)=\theta(x_{c1}-x)\frac{1}{2}x^{2}+\theta(x-x_{c1})\theta(x_{c2}-x)P(x)+\theta(x-x_{c2})T$
(182)

where $\theta(x)$ denotes the Heaviside step function, $P(x)$ denotes a third
order polynomial that ensures a smooth connection between the harmonic trap,
the potential barrier and the threshold potential $T$ whose role will be
discussed below. The coordinates $x_{c1}$ and $x_{c2}$ refer to connection
points, where the former connects the harmonic trap and the barrier and the
latter connects the barrier and the threshold. The bosons repel each other via
the contact interaction potential, and its strength is measured via the mean-
field parameter $\Lambda=\lambda_{0}(N-1)$.

Fig. 11 shows a schematic plot of the trapping potential for $t>0$ (solid
black). The maximum height of the barrier separates the IN and OUT regions,
denoted by the vertical red line. Moreover, the colored horizontal lines mark
the total energies of states with different numbers of bosons $N$ in the trap.
Depending on the threshold value $T$, certain states can become bound, meaning
that no additional particles can escape from the trap via tunneling through
the barrier. How the threshold can be utilized to control the tunneling-to-
open-space dynamics will be discussed below.

Figure 11: Schematic plot of the potential $V$ for $t>0$ (solid black). The
initial Gaussian-shaped density $\rho(x,t=0)$ is shown in solid blue. The IN
and OUT regions are separated at the position of the maximum of the potential
barrier (see solid red lines and labels). Colored horizontal lines denote the
energies of states with different numbers of bosons $N$ in the trap. The
energy differences between those states are denoted by the chemical potentials
$\mu_{i}$. The value $T$ denotes the height of the threshold potential in the
OUT region. All quantities are dimensionless. See text for details. The figure
is taken from Ref. [153].

Important quantities for the analysis are the density $\rho(x,t)$ [see Eq.
(2.4)], as well as the density in momentum space $\rho(k,t)$, or simply the
momentum distribution, which is obtained by performing a Fourier
transformation (FT) of the one-particle reduced density in real space.
Furthermore, the degree of fragmentation $f$, defined in Eq. (29), and
especially its evolution in time, is of major interest. The above quantities
are computed at the many-body level employing MCTDHB.

The analysis starts by looking onto the shape of the momentum distribution in
general. Fig. 12 shows $\rho(k,t)$ for a system of $N=101$ bosons with
interaction parameter $\Lambda=0.3$ at four different instants in time
$t_{1}<t_{2}<t_{3}<t_{4}$. It consists of two parts, a Gaussian background and
a peak structure. The former refers to the particles that did not escape from
the trap yet, since the initial real-space density is also a Gaussian. On the
contrary, the peaks refer to the momenta of the already escaped particles. The
heights of the peaks are time-dependent, which can be seen from the growth of
the dominant peak. The position of the peaks can be understood from a static
model introduced below.

Figure 12: The momentum distribution $\rho(k,t)$ for $N=101$ bosons with
repulsion strength $\Lambda=0.3$ at four instants in time with
$t_{1}<t_{2}<t_{3}<t_{4}$. The Gaussian background refers to the bosons
remaining in the trap, whereas the peaks refer to the momenta of the emitted
bosons. The height of the dominant peak is clearly growing in time. All
quantities are dimensionless. See text for details. The figure is adapted from
Ref. [152].

The minimal total energy $E_{\text{TOT}}$ of the system is comprised of the
energy of the bosons that remain in the IN region, denoted by
$E_{\text{HO}}(\lambda_{0},N_{\text{IN}})$ indicating its dependence on the
two-body repulsion strength $\lambda_{0}$ and the number of particles in the
trap $N_{\text{IN}}$, as well as on the minimal energy of the $N_{\text{OUT}}$
bosons that tunneled through the barrier into open space. Their minimal energy
is the height of the threshold $T$. One thus obtains

$E_{\text{TOT}}=E_{\text{HO}}(\lambda_{0},N_{\text{IN}})+N_{\text{OUT}}\,T.$
(183)

Assuming that the interaction among the emitted bosons is negligible, the
kinetic energy of the $i$-th escaped boson is given by

$E_{\text{kin}}(T,\mu_{i})=\mu_{i}-T$ (184)

where the $\\{\mu_{i}\\}$ are chemical potentials that refer to the energy to
bring an additional boson into the IN region with already $N-i$ bosons. The
corresponding momenta read

$k_{i}^{T}=\sqrt{2(\mu_{i}-T)}$ (185)

where $m=\hbar=1$ is assumed. It is stressed that this model describes the
tunneling of a single boson.

Figure 13: The tunneling to open space of $N=2$ bosons. Left panel: Static
consideration of the energy dependence of the possible final states
$|N_{\text{IN}},N_{\text{OUT}}\rangle$ on the threshold $T$. For $T\lesssim
0.5$, both particles are expected to tunnel, whereas for $0.5\lesssim
T\lesssim 0.8$, only one boson is expected to tunnel. For even larger $T$, the
bosons should remain trapped. Right panel: The obtained peaks in the momentum
distribution $\rho(k,t=600,T)$ for different $T$. For $T\leq 0.4$, two peaks
appear at $k_{1,2}^{T}$, for $T=0.5$ and $0.6$, only one peak appears. This is
in agreement with the expected number of ejected bosons from the trap. For
larger $T$, the peaks are shifted to smaller momenta since the bosons have to
overcome a larger threshold. All quantities are dimensionless. See text for
details. The figures are adapted from Ref. [153].

These model assumptions are tested against a system with $N=2$ bosons and two-
body repulsion strength $\lambda_{0}=1$. The left panel of Fig. 13 shows the
dependence of $E_{\text{TOT}}$ on the threshold value $T$ for the three
possible final states $|N_{\text{IN}},N_{\text{OUT}}\rangle$, given by
$|2,0\rangle$, $|1,1\rangle$ and $|0,2\rangle$. The utilized notation
$|N_{\text{IN}},N_{\text{OUT}}\rangle$ should not be confused with a Fock
state and is just used to count the particles in the two regions IN and OUT.
For $T\lesssim 0.5$, the lowest-in-energy state is $|0,2\rangle$, i.e., the
final state where both bosons have tunneled to open space. For $0.5\lesssim
T\lesssim 0.8$, the state $|1,1\rangle$ is the lowest, meaning that only one
boson is expected to tunnel in the long-time dynamics. For even higher $T$,
tunneling is suppressed completely because the state $|2,0\rangle$ corresponds
to the lowest energy. This is in full agreement with the results shown in the
right panel of Fig. 13. There, the peaks in the momentum distribution
$\rho(k,t=600)$, corresponding to the momenta $k_{1}^{T}$ and $k_{2}^{T}$ of
the emitted particles, are presented for different values of the threshold
$T$. The $k_{2}^{T}$-peaks are plotted upside-down for better visibility. The
first observation is that increasing $T$ shifts the peaks to smaller momenta,
which is due to the higher energy that is necessary to overcome the threshold.
Secondly, two peaks appear up to $T=0.4$, whereas for $T=0.5$ and $0.6$ only
one peak appears. This can be understood in terms of the static considerations
from the left panel, where it was observed that for $T\lesssim 0.5$ both
particles are expected to tunnel through the barrier, whereas only one boson
is expected to tunnel for $0.5\lesssim T\lesssim 0.8$. Thus, the threshold can
be utilized to control the number of tunneling bosons. The arrows in the
bottom of the right panel mark the calculated positions of the peaks with
respect to the model of Eq. (185). The good agreement with the positions of
the peaks of the momentum distribution suggests that the tunneling of single
particles can be well described by the above model.

Figure 14: Evolution of the natural occupations of the first few natural
orbitals. Left panel: The case of zero threshold, $T=0$, utilizing $M=4$
orbitals. The interaction strength is $\lambda_{0}=0.3$. For all particle
numbers shown, the system evolves from being condensed to being multi-fold
fragmented. Right panel: The occupation of the first two natural orbitals for
$N=2$ particles and $\lambda_{0}=1$, but for different heights of the
threshold $T$. In general, the higher the threshold the later the development
of fragmentation sets in. Moreover, the system with larger $T$ fragment less
than the ones with smaller $T$. All quantities are dimensionless. See text for
details. The figures are taken from Refs. [152, 153].

It remains the question whether the bosons do only tunnel one by one, or if it
is also possible that several bosons tunnel at once. To answer this question,
one has to consider that the peaks for $T\lesssim 0.5$ in Fig. 13 (right
panel) appear sequentially, i.e., first the dominant peak develops, and then
the second, smaller peak appears. The fact that the height of the peaks are
time-dependent suggests that the two processes, namely the separate ejection
of the first and second boson, are not independent. If however a multi-boson
ejection happened, peaks at higher momenta than the ones observed would
appear. For instance, in the limit of large $N$, the chemical potential
corresponding to a two-boson ejection would be $\mu^{2b}\approx 2\mu_{1}$,
which, for the case of $T=0$, yields $k^{2b,0}\approx\sqrt{2}k_{1}^{0}$. This
is simply not observed in the Fourier spectrum of $\rho(k,t)$. Thus, it can be
deduced that the bosons tunnel to open space one at a time. Nevertheless, as
mentioned above, the individual processes of the ejection of single bosons are
not independent, meaning that they quantum-mechanically interfere.

Finally, it is discussed whether the overall process of tunneling to open
space is indeed a many-body process or whether it can also be described at the
mean-field level. To this end, the time-evolution of the occupations of the
first few natural orbitals for systems with different numbers of bosons and
different heights of the threshold is studied in Fig 14. In the weakly-
interacting regime considered here, all initial states of the bosonic clouds
are essentially condensed. For the case of zero threshold (left panel), one
observes that for all numbers of bosons considered the system becomes multi-
fold fragmented. Especially for $N=4$ bosons, at least four orbitals become
macroscopically occupied. With respect to the case of $N=2$ particles and non-
zero thresholds, the behavior is similar, meaning that the initially coherent
BEC becomes highly fragmented. In general, the onset of fragmentation happens
at later times when $T$ is larger. Furthermore, the degree of fragmentation is
largest in the absence of any threshold.

To summarize, it has been found that the tunneling-to-open-space dynamics of
repulsive bosons from a harmonic trap are of many-body nature. With time, the
system evolves from being coherent to become highly fragmented, and thus a
mean-field description is inapplicable even in the weakly-interacting regime.
The bosons tunnel one by one, meaning that no multi-boson tunneling occurs
where more than a single boson at the same time is emitted jointly from the
trap. However, the individual tunneling processes are not independent, they
quantum-mechanically interfere. The threshold in the open space can be
utilized to control the number of emitted particles. In particular, it can be
used to suppress tunneling completely. A detailed analysis on the coherence of
the system, showing that the emitted bosons lose coherence with the bosons in
the trap and among themselves, can be found in Refs. [150, 151, 152, 153,
154], as well as additional details.

####  6.2.2 Phantom vortices

A recently found vortex type in a rotating and repulsive 2D single-component
BEC, termed phantom vortex, is presented in this section. Different to the
common understanding of a quantum vortex that, among other features, has a
density node (or core) and a phase discontinuity surrounding it, a phantom
vortex cannot be detected in the condensate density. Although they were not
found in single-component condensates so far, coreless vortices have been
predicted earlier for multi-component condensates, e.g., spinor condensates
where the core of one species, i.e., bosons with a certain spin, is filled by
another species [155, 156]. Since the experimental realization of gaseous BECs
[4, 6, 5, 157, 158], the nucleation of vortices in rotating BECs was of
interest and it was found that there is a critical rotation frequency below
which no quantized vortices appear [161, 159, 162, 163, 164, 160]. In contrast
to that, phantom vortices appear also below this critical rotation frequency.
Instead of a node in the density, they rather manifest themselves in nodes of
the densities of underlying natural orbitals. Thus, phantom vortices are a
many-body object that cannot be described at the mean-field level. Moreover,
their occurrence is accompanied by the onset of fragmentation. The following
results were published in Ref. [165]. Earlier results on phantom vortices in a
stirred BEC where, e.g., the fragmentation entropy is analyzed, can be found
in Ref. [166].

The system under consideration is a BEC of $N=100$ spinless bosons. They
interact via the same Gaussian repulsion as in Eq. (180) with interaction
parameter $\Lambda=\lambda_{0}(N-1)=17.1$. They are confined in a time-
dependent trap that is varied stepwise. It reads

$V(\mathbf{r},t)=\frac{1}{2}\left(x(t)^{2}+y(t)^{2}\right)+\frac{1}{2}\eta(t)\left(x(t)^{2}-y(t)^{2}\right)$
(186)

where $\eta$ is the parameter of the anisotropy given by the second term in
Eq. (186). The $x-$ and $y-$coordinates are varied in time according to

$\begin{pmatrix}x(t)\\\ y(t)\end{pmatrix}=\begin{pmatrix}\cos(\omega
t)&\sin(\omega t)\\\ -\sin(\omega t)&\cos(\omega t)\end{pmatrix}$ (187)

where the frequency $\omega$ is set to $0.78$. Initially, the trapping
potential is isotropic and harmonic, meaning that $\eta(t=0)=0$, and the
corresponding ground state is computed both at the mean-field level using the
GP equation as well as at the many-body level using imaginary time-propagation
of the MCTDHB equations. Afterwards, a slight anisotropy is added by
increasing $\eta$ from zero to its maximal value of $0.1$ until $t=80$. The
resulting elliptic trap is kept between $80<t\leq 300$, before $\eta$ is
ramped down to zero again between $300<t\leq 380$. For the last interval until
$t=500$, the isotropic harmonic trap is kept constant. The entire procedure is
shown in the upper panel of Fig. 15, together with representative densities of
the BEC for the individual intervals. These densities already suggest that at
no point in time there is a vortex contained in it, as discussed in more
detail below.

Figure 15: Top panel: Time-evolution of the anisotropy parameter $\eta(t)$
(red line) for the 4 distinct intervals discussed in the text, together with
representative densities (blue). The dynamic protocol is divided into a period
of ’ramp up’, ’maximal anisotropy’, ’ramp down’ and ’isotropic’. Second panel:
Time-evolution of the first 4 natural occupation numbers
$\rho_{i}^{(\text{NO})}(t)$. Fragmentation sets in at the end of the ramp up
of the anisotropy. Third panel: Time-evolution of the orbital angular momenta
$(L_{z})_{ii}(t)$ and the total angular momentum per particle, $L_{z}(t)/N$.
Fluctuations are maximal when fragmentation sets in. Bottom panel: Comparison
between the time-evolution of the total energy $E$ relative to the initial
energy $E_{0}$ and $L_{z}(t)/N$. The GP theory underestimates both quantities,
especially the angular momentum is approximately three times less than at the
many-body level. All quantities are dimensionless. See text for details. The
figure is taken from Ref. [165].

Of particular interest for the analysis will be, apart from the density
$\rho(\mathbf{r})$, the natural occupation numbers, denoted here by
$\rho_{i}^{(\text{NO})}$ but with the same meaning as the $n_{i}$ in Eq. (28),
and the densities $|\phi_{i}|^{2}$ of the natural orbitals (mind the different
notation as compared to Section 3 where the $\phi_{i}$ denote the working
orbitals of MCTDHB). Furthermore, the time-evolution of the expectation values
of the orbital angular momenta, given by
$(L_{z})_{ii}=\langle\phi_{i}|\hat{L}_{z}|\phi_{i}\rangle$, and of the total
angular momentum per particle,
$L_{z}/N=\frac{1}{N}\sum_{i,j=1}^{M}\rho_{ij}\,(L_{z})_{ij}$, are analyzed.
The time-argument is suppressed in the above mentioned quantities.
Computations are carried out for orbital numbers $M=4$ (many-body) and $M=1$
(GP). Additionally, signatures of phantom vortices can also be seen in the
phase of the first-order correlation function
$g^{(1)}(\mathbf{r}|\mathbf{r}^{\prime};t)$ [see Eq. (D.6) in Appendix D]
defined by

$S_{g}(\mathbf{r}|\mathbf{r}^{\prime};t)=\text{arg}[g^{(1)}(\mathbf{r}|\mathbf{r}^{\prime};t)]$
(188)

as well as in the phases of the natural orbitals given by

$S_{i}(\mathbf{r},t)=\text{arg}[\phi_{i}(\mathbf{r};t)].$ (189)

The discussion of the results is split into two parts. At first, the
appearance of phantom vortices is analyzed at an instant in time when the trap
is already isotropic again. Afterwards, the degree of fragmentation in the
system is discussed.

Figure 16: Top panels: Total density $\rho(\mathbf{r})$ (a) and densities of
the natural orbitals (b)-(e) at $t=450$, i.e., when the trap is already
isotropic again. Although no vortex is observed in the total density, phantom
vortices in the orbitals $\phi_{1}$, $\phi_{3}$ and $\phi_{4}$ appear. Bottom
panels: Phases of the first-order correlation function (f) and of the natural
orbitals (g)-(j). Phase discontinuities are observed for all phantom vortices.
All quantities are dimensionless. See text for details. The figure is taken
from Ref. [165].

With respect to the first part, the total density and the densities of the
natural orbitals at $t=450$, i.e., when the anisotropy is ramped down to zero
again, are examined in Fig. 16. One observes that the total density [panel
(a)] does not exhibit any vortex. On the contrary, several phantom vortices
appear in the densities of the natural orbitals. The densities of the first
and fourth natural orbitals contain a charge-1 and a charge-3 phantom vortex,
respectively [panels (b) and (e)]. The density of the third orbital contains a
pair of charge-1 phantom vortices [panel (d)]. The mechanism responsible for
the nucleation of the phantom vortices in $\phi_{1}$ and $\phi_{4}$ is
different to the one responsible for the pair in $\phi_{3}$. The former
mechanism is called node mutation, and is explained for the charge-1 phantom
vortex in $\phi_{1}$. This phantom vortex originates from the node of the
initial second orbital $\phi_{2}$, which resembles a 2$p_{x}$ orbital that has
an angular node. When angular momentum is added to the system due to the
anisotropy, the lobes of the $2p_{x}$ orbital spread out and tend to close
this angular node, leaving an elliptic density node in the center. After the
trap is isotropic again, a phantom vortex is left at the center, which now
appears in the density of $\phi_{1}$ since the first and second orbital have
switched labels at $t\approx 220$ because of their changing occupation numbers
[see second panel from the top in Fig. 15]. Thus, the remaining phantom vortex
was nucleated from an initial node. A similar process happens for the charge-3
phantom vortex in the density of $\phi_{4}$. It originates from
$\phi_{3}(t=0)$ which resembles a $2p_{y}$ orbital. There, the angular node
mutated first to three single phantom vortices that finally merged into a
larger charge-3 phantom vortex. The fact that it finally appears in the fourth
natural orbital is due to the label switching of $\phi_{3}$ and $\phi_{4}$
during the ramp-up of the anisotropy.

Figure 17: Time-evolution of the fragmentation for different numbers of
particles $N$. The same dynamical protocol of the trapping potential as in
Fig. 15 is used. For all considered boson numbers, the occupation
$\rho_{1}^{(\text{NO})}$ decreases clearly. In general, fragmentation sets in
later for higher particle numbers. Results are computed using $M=2$ orbitals.
The interaction parameter is kept fixed at $\Lambda=17.1$. All quantities are
dimensionless. See text for details. The figure is taken from the
supplementary material of Ref. [165].

The mechanism responsible for the phantom-vortex pair in $\phi_{3}$ [panel
(d)] is called slow orbital-orbital vorticity transfer because of the
following observation. At $t\approx 220$, the pair of co-rotating phantom
vortices appeared in the density of $\phi_{2}$. When time progressed, the two
phantom vortices started to move towards the density edge and disappeared
finally at $t\approx 275$. At the same time, two phantom vortices were build
at the density edge of $\phi_{3}$, which afterwards continued to move towards
the center. At $t\approx 340$, the separation between this pair is the same as
before in the density of $\phi_{2}$. Thus, vorticity has been transferred
between the two natural orbitals over a period much longer than the rotation
period of $\tau=2\pi/\omega\approx 8$. This process between the second and
third orbital can also be seen from the correlated orbital angular momenta
$(L_{z})_{22}$ and $(L_{z})_{33}$ in Fig. 15 for $t\geq 250$. The position of
the phase discontinuities shown in Figs. 16(g), (i) and (j) coincide with the
corresponding positions of the phantom vortex cores. Especially in panel (j)
one can see from the merging three discontinuity branches in the center that
the corresponding phantom vortex in panel (e) is of charge 3. Comparing the
density in panel (a) and the phase $S_{g}(\mathbf{r}|0;450)$ in panel (f), one
observes the occurrence of several so-called ghost vortices at the density
edge. The latter are distinct from phantom vortices because they do not enter
the bulk density of the condensate and do also not contribute much to the
total angular momentum.

The dynamical fragmentation that occurs during the time-evolution, as well as
its consequences for the angular momenta and the energy of the system [see
again Fig. 15] is now analyzed. The initial state is essentially condensed
($\rho_{1}^{(\text{NO})}\approx 99.7\%$), and it stays condensed for almost
the entire period in which the anisotropy is ramped up to its maximal value.
In this time frame, the total angular momentum and energy computed at the
many-body and mean-field levels coincide as long as the BEC remains coherent.
Once fragmentation sets in, the latter quantities start to deviate from each
other. During the first part of the time interval in which the trap is kept
maximally anisotropic ($80<t\leq 150$), $\rho_{1}^{(\text{NO})}$ decreases
significantly to roughly $40\%$, and the second, third and fourth natural
orbitals become macroscopically occupied. At the same time, the orbital
angular momenta vary the most and several of them reach their maximum values.
Then, the occupations and angular momenta of $\phi_{1}$ and $\phi_{4}$ remain
more or less constant. For the orbitals $\phi_{2}$ and $\phi_{3}$, both
quantities change slightly in a correlated manner, i.e., when
$\rho_{2}^{(\text{NO})}$ increases $\rho_{3}^{(\text{NO})}$ decreases by the
same amount. This again reflects the slow orbital-orbital vorticity transfer
between $\phi_{2}$ and $\phi_{3}$ discussed above. With regard to the total
angular momentum and energy, both quantities vary strongly during the periods
of maximal anisotropy and ramp down, and saturate only afterwards. The
corresponding GP values are substantially different since they underestimate
the system’s gain of energy and angular momentum. In particular for the latter
quantity, it stays below $0.5$, well below the threshold for the occurrence of
vortices at the mean-field level at unit angular momentum per particle.
Interestingly, $L_{z}/N$ for $M=4$ orbitals is greater than this threshold,
but still no vortices are observed in the density $\rho(\mathbf{r})$. This can
be rationalized because one can show analytically that a vortex in the density
of a fragmented BEC, i.e., a fragmented vortex, requires more than unit
angular momentum per particle [see the supplementary information of Ref.
[165]]. The reason is that a fragmented vortex originates from coincident
phantom vortices, meaning that all natural orbitals must have a phantom vortex
at the same position, and all of these phantom vortices need to have different
charges such that the natural orbitals are orthogonal. As seen from panels
(b)-(d) in Fig. 16, there are no coincident phantom vortices, and thus the
fragmented system does not show any vortices in the density.

The onset of fragmentation was also observed for different numbers of bosons.
Fig. 17 shows the evolution of $\rho_{1}^{(\text{NO})}$ for different $N$
between $10$ and $10^{4}$, utilizing $M=2$ time-adaptive orbitals. The
interaction parameter $\Lambda=17.1$ is kept fixed. One generally observes
that fragmentation sets in for all considered particle numbers, and that
furthermore the onset of fragmentation occurs at later times for higher
particle numbers.

In conclusion, the nucleation of phantom vortices has been investigated for a
rotating BEC that acquires angular momentum via a time-dependent modification
of the trap. Two different mechanisms for the creation of phantom vortices
were found, namely node mutation as well as slow orbital-orbital vorticity
transfer. Phantom vortices, by definition many-body objects, cannot be
described at the mean-field level. The initially coherent system highly
fragments, also for different numbers of particles in the condensate. Further
information can be found in Ref. [165], in particular in the supplementary
material containing videos that show, e.g., the evolution of all quantities
shown in Fig. 16.

###  6.3 Excitation spectra

In this section, details on current applications of LR-MCTDHB are presented.
The amount of scientific publications where it is utilized is still rather
small because the theory itself has been developed only a few years ago [34,
35]. In addition to that, the efficient MPI-parallelized numerical
implementation of LR-MCTDHB was a very demanding task, and all the technical
challenges going along with it, extensively described in Section 4, were
solved very recently. It is stressed again that without such an
implementation, a vast majority of problems cannot be treated at all,
especially in $D>1$ spatial dimensions. Before developing the sophisticated
and complex parallel implementation, there was a first implementation
available that could treat small 1D systems. The results of Refs. [34] and
[167] were obtained with this code [see below and Section 6.3.1]. In more
recent applications, discussed in Secs. 6.3.2 and 6.3.3, the new
implementation was utilized since it was impossible to obtain the same results
with the former code.

(a)

(b)

Figure 18: (a) Numerical convergence of the first three center-of-mass excited
states (labeled by 1, $2^{\prime}$ and $3^{\prime}$) with respect to the
number of orbitals $M$ for $N=10$ bosons in a harmonic trap with trap
frequency $\omega_{\text{ho}}=\sqrt{2}$. The strength of the contact repulsion
is $\lambda_{0}N=1$. For $M=5$ orbitals, the energies obtained from LR-MCTDHB
coincide with the analytic predictions. The numerical implementation used is
the precursor of the implementation described and benchmarked in Sections 4
and 5. (b) Comparison between the excitation spectra obtained from BdG
(bottom) and LR-MCTDHB($2$) (top) for $N=10$ bosons in a shallow double well.
The interaction strength is the same as in panel (a). Excitations labeled with
primes refer to states solely obtained at the many-body level. The intensity
of the response to a linear perturbation ($f^{+}=f^{-}=x$) is denoted by red
triangles, whereas the intensity of the response to a quadratic perturbation
($f^{+}=f^{-}=x^{2}$) is denoted by green squares. The insets show the trap
potential and the ground-state densities. One observes that the many-body low-
energy spectrum contains a lot more states than the mean-field spectrum. All
quantities are dimensionless. See the main text and Ref. [34], from which both
figures are taken, for further details, especially for the discussion of the
excitation type shown at the top of panel (b).

The main purpose of Ref. [34] is to introduce LR-MCTDHB by giving a derivation
of its main equations, presented in Section 3.2.2 of this work. In particular,
the equations for the submatrices as described in Eq. (93) are developed for
the case of contact interaction between the bosons. As an application,
excitations for BECs in harmonic and shallow double-well traps are analyzed.
Fig. 18(a) shows the numerical convergence for the first three c.m. excited
states (labeled by 1, $2^{\prime}$ and $3^{\prime}$) of $N=10$ bosons in a
harmonic trap with trap frequency $\omega_{\text{ho}}=\sqrt{2}$ in terms of
the number of orbitals $M$. The strength of the repulsion is given by
$\lambda_{0}N=1$. One can observe that the c.m. excitations are accurately
obtained for $M\geq 5$ orbitals. The BdG theory only yields the first c.m.
excitation. In Fig. 18(b), the low-energy spectrum as well as the response
intensities for linear and quadratic perturbations ($f^{+}=f^{-}=x,\,x^{2}$)
are presented for $N=10$ bosons in a shallow double well of the form
$V(x)=b/2\cos\left(\pi x/3\right)+\omega_{\text{ho}}^{2}x^{2}/2$ with barrier
height $b=5$. The most important observation by comparing the number of
obtained states for BdG and LR-MCTDHB($2$) is that the many-body spectrum
contains a lot more states. This in turn means that although the degree of
depletion is only about $0.2\%$, a many-body description for the lowest-in-
energy excitation spectrum is unavoidable in order to obtain accurate results.
In the subsequent Section 6.3.1, the excitation spectrum for the case of BECs
trapped in harmonic and double-well potentials is elaborated in more detail,
especially with respect to the question of how many-body excitations can be
triggered in the out-of-equilibrium dynamics of the condensates.

Before discussing further applications of LR-MCTDHB in the following, a remark
on how mean-field and many-body excitations are defined in this work is made
in order to avoid misconceptions of these important terms. Excitations
calculated from Eq. (116), i.e., from LR-MCTDHB ($M>1$), are termed many-body
excitations since they stem from a many-body theory. On the contrary, all
excitations calculated with the BdG equations in Eq. (A.16) of Appendix A.1.2
are referred to as mean-field excitations since they stem from a mean-field
theory. Moreover, excitations where only a single boson is excited from the
condensed mode are called single-particle excitation, whereas excitations
where more than a single boson is excited are called multi-particle
excitations. These definitions are utilized consistently in the applications
below.

####  6.3.1 One-dimensional harmonic and double-well systems

As a first application of LR-MCTDHB, the many-body nature of excitations in a
1D harmonic trap as well as in shallow and deep symmetric double wells is
analyzed. The low-energy spectra obtained from the BdG equation at the mean-
field level and from LR-MCTDHB at the many-body level are compared.
Furthermore, it is investigated which excitations are involved in the dynamics
due to certain quench scenarios. The following results where published in Ref.
[167].

Excited states are studied in a harmonic trapping potential with a Gaussian
barrier given by

$V(x)=ax^{2}+b\,\exp\left(-cx^{2}\right)$ (190)

where $a$ denotes the harmonic trap frequency and $B$ denotes the barrier
height. For simplicity, $\hbar=m=1$ is assumed. Furthermore, the trapping
parameters are set to $a=1/2$ and $c=1$. The bosons interact via the zero-
ranged contact interaction potential. The interaction strength is again
measured by the mean-field parameter $\Lambda=\lambda_{0}(N-1)$ where
$\lambda_{0}$ is chosen to be positive to account for repulsion. Fig. 19 shows
$V(x)$ for different barrier heights, together with the corresponding ground-
state densities. For the purely harmonic trap ($b=0$), the BEC is essentially
condensed with $n_{1}=99.9\%$. On the contrary, for the shallow ($b=5$) and
deep ($b=10$) double wells, the ground state is already fragmented by
approximately $5\%$ and $40\%$, respectively. One can further observe that for
the deep well, the overlap between the densities in the left and right wells
is marginal, whereas there is still a clear overlap observable in the shallow
case.

Since the trap is reflection invariant, the excitations can be categorized in
gerade ($g$) and ungerade ($u$) symmetry. Fig. 20 shows the low-energy
spectrum of $N=10$ bosons for different barrier heights $0\leq b\leq 10$ and
$\Lambda=1.0$. Excitation energies are given relative to the ground-state
energy, i.e., $\Delta E=E-E_{0}$. One readily observes the formation of bands,
which is even more pronounced for lower interaction strengths [see Fig. 2 in
Ref. [167]]. The bands are labeled by integer numbers, where the zeroth band
denotes the one that is lowest in energy. Comparing the mean-field and many-
body approaches, one can see that the BdG spectrum only contains 5 different
excitations, whereas the spectrum computed with LR-MCTDHB shows a much richer
structure where a lot more excited states appear. In terms of numerical
convergence, $M=2$ orbitals are sufficient to describe the spectrum up to
$\Delta E\sim 2$ for not too low barrier heights, i.e., $b\geq 4$. For states
higher in energy, at least four orbitals are necessary.

Figure 19: Trap potentials (solid red) and ground-state densities (gray) for
different barrier heights $b$. The ground states are computed with
MCTDHB($2$). The BEC consists of $N=10$ bosons and the repulsion strength is
$\Lambda=1.0$. The occupation number $n_{1}$ of the first natural orbital
shows that increasing the barrier height leads to fragmentation. Vertical
lines refer to the lowest-in-energy excited states computed with LR-
MCTDHB($2$). The green potentials and the blue eigenvalues refer to the
harmonic approximation of the right well. All quantities are dimensionless.
See text for details. The figure is taken from Ref. [167]. Figure 20: Low-
energy excitation spectrum of $N=10$ bosons with interaction parameter
$\Lambda=1.0$ for different barrier heights $b$. Many-body results computed
with LR-MCTDHB($2$) are denoted by solid lines (gerade excitations in green,
ungerade excitations in red). Results for $M=4$ orbitals are given by dotted
lines, and the corresponding BdG energies are given by open black circles. All
states up to $\Delta E\sim 2$ are converged with two orbitals. The BdG
approach misses many states in the shown energy range. All quantities are
dimensionless. See text for details. The figure is taken from Ref. [167].

In the following, different shift and quench protocols are applied to trigger
excitations from the low-energy spectrum. For gerade excitations, a quench of
the harmonic confinement from $a=0.4\rightarrow 0.5$ is performed, whereas
ungerade excitations are addressed by a shift of the trap position, meaning
that $x\rightarrow x+x_{\text{shift}}$ with $x_{\text{shift}}=0.1$. The gerade
excitations are referred to as breathing modes, and the ungerade excitations
are referred to as dipole modes. The dynamics of the bosonic clouds for the
two cases are shown in Fig. 21 on a short time scale. In the left panel,
dipole oscillations for a BEC in the deep double well are depicted, which
means that in each well, the bosons oscillate around the equilibrium position
given by the minima of the wells. On the RHS, the breathing dynamics for
$N=10$ bosons in the purely harmonic trap is shown. One can observe that the
density spreads to the left and right in a symmetric manner, and returns to
the minimum of the well. The Fourier spectra of the time-evolution of certain
quantities is used to extract the spectrum of involved excitations in the
dynamics. This comprises $(a)$ the expectation value $\langle x\rangle$, $(b)$
the variance per particle $\frac{1}{N}\Delta^{2}_{\hat{X}}$, and $(c)$ the
density $\rho$ at the position $x=1$. The many-body variance is discussed in
detail in Appendix D. In Ref. [167], the Fourier spectra for $b=0\,,5$ and
$10$ are investigated. Here, the discussion is restricted to the latter. To
remind the reader, the system is highly fragmented by approximately $40\%$ in
this case.

Figure 21: Many-body dynamics of the density of $N=10$ bosons with interaction
parameter $\Lambda=1.0$ for a shift of the trap position in the deep double
well (left panel) and a quench of the harmonic confinement (right panel). In
the former scenario, one observes dipole oscillations of the bosonic clouds in
the two wells, whereas one observes breathing oscillations for the dynamics
due to the quench of the confinement. All quantities are dimensionless. See
text for details. The figure is taken from Ref. [167].

Table 4 shows excitation and de-excitation energies where the obtained values
from the BdG and LR-MCTDHB($2$) theories (columns 2 and 3) are compared to the
obtained excitation energies of the FT spectra from the shift and quench
scenarios (columns 4-7). Moreover, several de-excitations, i.e., transitions
between a higher and a lower excited state, are presented (last 4 rows). Those
refer to higher-order effects and are in principle not obtained directly by
LR. However, since LR yields the exact excitation spectrum once it is applied
to an exact eigenstate of the full many-body Hamiltonian, the de-excitations
can be extracted from differences between the excitation energies of the LR
spectrum. The excitation labels are shown in the first column. Many-body
excitations, i.e., excitation that only appear in the many-body spectrum, are
denoted by primed labels.

|  $M=1$ |  $M=2$ |  Shift GP |  Shift MB |  Quench GP |  Quench MB
---|---|---|---|---|---|---
$0u$ |  0.035 (0.46) |  0.063 (0.19) |  0.034 (0.02a) |  |  |
$0g^{\prime}$ |  |  0.064 |  |  |  |  0.066 (0.36b)
$1g$ |  2.085 (0.91) |  2.092 (0.81) |  |  |  2.085 (1.0b,0.93c) |  2.092 (0.33b,0.89c)
$1u$ |  2.102 (0.46) |  2.095 (0.57) |  2.102 (0.94a,0.66c) |  2.095 (0.93a,0.7c) |  |
$1u^{\prime}$ |  |  2.158 (0.08) |  |  2.163 (0.02a,0.02c) |  |
$1g^{\prime}$ |  |  2.161 (0.09) |  |  |  |  2.161 (0.12b)
$2g$ |  3.912 (0.04) |  3.934 (0.03) |  3.914 (0.04c) |  |  3.914 (0.07c) |  3.934 (0.07b,0.1c)
$2u$ |  3.997 (0.07) |  3.943 (0.07) |  3.997 (0.04a,0.23c) |  3.943 (0.03a,0.15c) |  |
$2u^{\prime}$ |  |  4.029 (0.05) |  |  4.032 (0.02a,0.07c) |  |
$2g^{\prime}$ |  |  4.039 |  |  |  |
$1g^{\prime}\leftrightarrow 1g$ |  |  0.069 |  |  |  |
$2g\leftrightarrow 1g$ |  1.827 |  1.841 |  |  |  |  1.893 (0.01b)
$2u\leftrightarrow 1u$ |  1.894 |  1.848 |  1.895 (0.01c) |  1.849 (0.02c) |  |
$1g\leftrightarrow 0g^{\prime}$ |  |  2.028 |  |  |  |  2.032 (0.02b)
Table 4: Table of excitation and de-excitation energies for $N=10$ bosons with
$\Lambda=1.0$ in the deep double well. Listed are results from BdG and LR-
MCTDHB($2$), as well as from the Fourier spectra of the expectation value of
the position $\langle x\rangle$ (superscript $a$), the variance per particle
$\frac{1}{N}\Delta^{2}_{\hat{X}}$ (superscript $b$) and the real-space density
$\rho(x=1)$ (superscript $c$), for both the shift and quench scenarios. GP
denotes the mean-field and MB the many-body case. The intensities are given in
brackets. All quantities are dimensionless. See text for more details. The
table is adapted from Ref. [167].

Up to $\Delta E\sim 4$, the many-body approach yields ten excited states,
whereas on the contrary the BdG spectrum only consists of five excitations.
The energies of states that appear in both spectra differ slightly from each
other for $1g$, $1u$, $2g$ and $2u$, but significantly for the lowest
excitation $0u$. That demonstrates that strong deviations between the mean-
field and many-body excitation energies can already occur for the lowest
excited states. Also the intensities of several states differ clearly from
each other, especially for the lowest states $0u$, $1g$ and $1u$. The
intensities from the LR calculations are identical to the response weights
described in Eq. (3.2.2), with $f^{+}=f^{-}=x$ for the ungerade and
$f^{+}=f^{-}=x^{2}$ for the gerade states.

Figure 22: Fourier spectra of $\langle x\rangle$,
$\frac{1}{N}\Delta^{2}_{\hat{X}}$, and $\rho(x=1)$ from the shift and quench
dynamics in the deep double well. The GP spectra are denoted by solid black
and the MCTDHB($2$) spectra by dashed red lines. The intensities are
normalized by the sum of the intensities of all constituent peaks, and
excitations with less than $1\%$ intensity are omitted. Especially the
spectrum of the many-body variance consists of several peaks that are not
observed in the corresponding GP spectrum. All quantities are dimensionless.
See text for details. The figure is adapted from Ref. [167].

Columns 4 and 5 show the obtained peaks in the Fourier analysis of $\langle
x\rangle$ and $\rho(x=1)$ for the shift scenario, both for the GP and many-
body cases. The Fourier spectra are shown in the left panels of Fig. 22.
Surprisingly, the shift also triggered a gerade excitation ($2g$) at the mean-
field level, whereas at the many-body level only ungerade states are observed.
Moreover, no peak is observed for $0u$ in the many-body dynamics, presumably
because its intensity is too weak to be distinguished from the background
noise. The obtained peak positions coincide with the predictions from BdG and
LR-MCTDHB($2$), respectively. With regard to de-excitations, the mean-field
and many-body approaches yield a very small peak for the $2u\leftrightarrow
1u$ transition in the density spectrum.

The Fourier analysis of the quench dynamics of
$\frac{1}{N}\Delta^{2}_{\hat{X}}$ and $\rho(x=1)$ are shown in column 6 (GP)
and column 7 [MCTDHB($2$)]. Only peaks corresponding to gerade excitations are
observed. Both approaches yield peaks for $1g$ and $2g$. At the many-body
level, two additional peaks referring to $0g^{\prime}$ and $1g^{\prime}$ are
obtained. The peak positions again accurately match the predictions from the
static analysis. Most importantly, the intensities of these excitations in the
variance spectrum are comparatively large, and therefore $0g^{\prime}$ and
$1g^{\prime}$ are the most promising many-body excited states to be detected
in an experiment. Furthermore, no de-excitation peaks are found in the GP
spectra, while the variance spectrum shows small peaks for the
$2g\leftrightarrow 1g$ and $1g\leftrightarrow 0g^{\prime}$ transitions at the
many-body level, coinciding with the values from LR-MCTDHB($2$).

In general, the sensitivity of the utilized operators to the detection of
excited states can be seen as follows. The spectrum of the local operator
$\rho(x=1)$ is in particular useful to detect ungerade excitations in the
shift scenario. The intensities of the peaks of many-body excitations are
however weak. To detect gerade excitation, nonlocal quantities, especially the
variance, turn out to be more sensitive than the density. In addition, the
variance spectrum yields intense peaks also for many-body excitations.

To summarize, it has been demonstrated that from a methodological point of
view, utilizing LR-MCTDHB in order to obtain the many-body excitation spectra
of repulsive BECs in 1D harmonic and double-well trap geometries is a
promising and highly accurate technique. Its predictions coincide with the
Fourier spectra of the time-evolution of different quantities like the local
density or the position variance. From a physical point of view, some of the
obtained excited states corresponding to many-body states not included in
standard mean-field approaches yield intense responses especially in the
variance spectra and are therefore likely to be measured in experiments. More
details on the presented results can be obtained from Ref. [167].

####  6.3.2 Triple wells and larger lattices

In this section, many-body effects in the excitation spectra of weakly-
interacting BECs in 1D lattices are investigated. The main objective is to
demonstrate with numerically converged results that a marginal depletion of
the order of $1\%$ is sufficient in order to observe clear differences between
the lowest-in-energy excitations computed at the mean-field and many-body
levels. The results presented below were recently published in Ref. [168].

Studying dynamics and excitations of ultracold bosons in 1D optical lattices
was of high interest in the past, and the literature is extensive in this
research field. Examples are the dynamics in the superfluid phase [169],
nonlinear dynamics in periodic potentials [170] or in fluctuation-driven
binary condensates [171]. A review of earlier works can be found in Ref.
[172]. Of particular relevance for the findings discussed in this section are
studies were the BdG theory does not agree with numerical and experimental
observations. It has been found that the oscillation frequency of a superfluid
BEC in a large optical lattice with a superimposed harmonic confinement
deviates clearly from the BdG prediction [173, 174]. In addition, the strong
depletion of a gaseous BEC in optical lattices cannot be fully explained by
the BdG theory, not only in 1D [175].

Optical lattices, in particular triple-well systems, have also been analyzed
recently utilizing MCTDHB. This includes the out-of-equilibrium dynamics of
rather small BECs consisting of a few particles only, induced by quenching the
interaction and the lattice depth [176, 177, 178] and periodic driving of the
system [179, 180]. Of particular interest was the analysis of cradle modes,
dipole oscillations corresponding to the transport of bosons over the barrier
between adjacent wells. Furthermore, the many-body dynamics of BECs initially
prepared in Mott-insulating states have been of interest [181]. The different
phases of the BEC, i.e., the superfluid, the Mott-insulating, and the
fermionized gas phase, can be detected from the many-body entropy and the
Glauber correlation functions [182]. It was even shown for a BEC in a tilted
triple well with hard-wall boundary conditions how the correlations between
wells can be controlled by the well depth, the interaction strength and the
tilt [183]. Many-body effects in the ground state of a bosonic mixture in a
more complex 2D lattice where studied using the multi-layer version of MCTDHB
[184].

Many-body excitations in optical lattices were addressed as well, either by
making a Fourier analysis of certain quantities [176] or by doing an exact
diagonalization. There, the system’s wave function was expressed using a many-
body ansatz based on a number state expansion including localized Wannier
functions [179, 180, 177, 178]. The latter approach of exact diagonalization
is possible for small systems with a few bosons only.

(a)

(b)

(c)

(d)

Figure 23: Low-energy spectra of $N=10$ bosons in a triple well computed with
LR-MCTDHB($7$). (a) Shallow case with lattice depth $V_{0}=1.01\,E_{R}$. For
$\Lambda=0$, the excitations form levels of degenerate states. Degeneracies
are lifted for $\Lambda>0$. The BdG results only contain the two doublets
denoted by BdG(1) and BdG(2). (b) Same as in (a) but in the deep lattice with
$V_{0}=10.13\,E_{R}$. The BEC for $\Lambda=2.0$ is fragmented by $8.4\%$. (c)
Comparison between BdG(1) and its multiples (open circles) and the many-body
results (colored closed symbols) for the system in (a) with $\Lambda=4.0$
where only the states $(+1)$ and $(-1)$ of the first single-particle band are
occupied by $N^{(1)}$ bosons. The difference grows with $N^{(1)}$. Inset:
Evolution of the relative error as defined in the main text. (d) Same as in
(c), but for the system in (b) with $\Lambda=2.0$. All quantities are
dimensionless. See text for details. The figures are adapted from Ref. [168].

Below, the lowest-in-energy excitations of BECs in 1D lattice potentials with
periodic boundary conditions are investigated using LR-MCTDHB. This opens the
possibility to study also larger systems with more wells and particles. The
utilized trap potential reads

$V(x)=V_{0}\,\cos^{2}\left(\frac{\pi}{l}\,x\right)$ (191)

where $V_{0}$ denotes the lattice depth and $l$ denotes the distance between
two neighboring lattice sites. The depth of the lattice is commonly expressed
in terms of the recoil energy $E_{R}=\frac{\hbar^{2}k_{0}^{2}}{2m}$ with the
lattice momentum $k_{0}=\frac{\pi}{l}$. For simplicity, $\hbar=m=l=1$ is
assumed. The bosons repel each other via the contact interaction potential and
its strength is again expressed in terms of the mean-field parameter
$\Lambda=\lambda_{0}(N-1)$.

(a)

(b)

Figure 24: (a) Same as in Fig. 23(a) but in a lattice with 10 sites. The
depletion grows faster with $\Lambda$ than in the triple well, meaning that
the system is even more sensitive for many-body effects. (b) Same as in Fig.
23(c) but for 10 sites and $\Lambda=1.0$. The relative error (inset) reaches
similar values as for systems in the shallow triple well with comparable
depletion. All quantities are dimensionless. See text for details. The figures
are adapted from [168].

As a start, $N=10$ bosons confined in a triple-well potential with periodic
boundary conditions are considered. The low-energy excitation spectra for both
a shallow ($V_{0}=5.0\approx 1.01\,E_{R}$) and a deep well ($V_{0}=50.0\approx
10.13\,E_{R}$) are discussed. Results are presented in Figs. 23(a) and (b). In
both wells, the spectra for the non-interacting BEC are structured in
sequences of levels with different degeneracies. The states and levels can be
associated with appropriate quantum numbers, envisioning that for $\Lambda=0$
eigenstates of a single particle in a periodic lattice are exact quasimomentum
eigenstates. The eigenstates appear in a band structure that is more
pronounced for deep lattices. For the shallow lattice in panel (a), the first
level consisting of two degenerate states denotes the excitations of putting
one boson in either a state with quasimomentum $(+1)$ or with quasimomentum
$(-1)$. Similarly, the states of level 5 denote the excitations where one
boson is put to either $(+1)$ or $(-1)$ from the second single-particle band.
To this end, the notation
$\left(n_{+1}^{(1)},n_{-1}^{(1)};n_{+1}^{(2)},n_{-1}^{(2)}\right)$ is
introduced, where the first two entries in brackets denote the occupation of
$(+1)$ and $(-1)$ from the first single-particle band and the last two entries
denote the occupation of $(+1)$ and $(-1)$ from the second single-particle
band. As another example, level 2 contains the states $(2,0;0,0)$, $(0,2;0,0)$
and $(1,1;0,0)$. It is stressed that the BdG spectrum only contains the
excitations of levels 1 and 5, and misses all multi-particle excitations in
the low-energy part. However, for $\Lambda=0$, the missing levels can be
anticipated by taking sums and multiples of the energies of the first and
fifth levels. As soon as there are correlations between the bosons, i.e., for
$\Lambda>0$, this is no longer possible, as discussed below. In panel (b),
only the first single-particle band appears in the shown energy range because,
as shown in Fig. 25, the first and second band are strongly separated from
each other.

Figure 25: Proof of numerical convergence for the systems discussed in Figs.
23 and 24. LR-MCTDHB($7$) clearly suffices in all cases. For the BEC in the
deep triple well [middle panel], all 65 single- and multi-particle excitations
below the bottom of the second band [denoted by BdG(2)] are obtained. Already
for $M=3$, the results are highly accurate which is due to the large
separation of the first two single-particle bands. All quantities are
dimensionless. See text for details. The figures are adapted from Ref. [168].

For non-zero repulsion in the shallow well in panel (a), the BECs are slightly
depleted ($f=1.1\%$ for $\Lambda=4.0$), whereas for the deep well in panel (b)
the BEC for $\Lambda=2.0$ is already fragmented by $8.4\%$. One can observe
that all excited states grow in energy once the bosons interact, and that
splittings of several states from the same level occur although they were
degenerate for $\Lambda=0$. Thus, the clear level structure of the non-
interacting case vanishes as $\Lambda$ increases, which makes it more
difficult to identify the individual states with regard to the categorization
in terms of quasimomentum states. The corresponding mean-field results contain
only two doublets for the shallow and one doublet for the deep well, denoted
by BdG(1) and BdG(2). All other states, representing multi-particle
excitations where more than one boson at a time is excited from the ground
state, are absent.

Figs. 23(c) and (d) demonstrate in how far the BdG energies can be utilized to
anticipate the positions of the missing multi-particle states. As mentioned
above, this yields the exact energies for the missing states in the non-
interacting system. In particular, multi-particle excitations of putting in
total $N^{(1)}$ bosons into the states $(+1)$ and $(-1)$ of the first single-
particle band are discussed. With regard to the BEC with $\Lambda=4.0$ in the
shallow triple well [panel (c)], one observes that BdG overestimates the
excitation energies for the states of $1\leq N^{(1)}\leq 5$, i.e., the
accurate many-body values are lower. The discrepancy between mean-field and
many-body results grows with $N^{(1)}$. Especially the relative error, defined
as
$E_{\text{rel}}=\left|\frac{\omega_{\text{BdG}}\left(N^{(1)}\right)-\omega_{\text{MB}}\left(N^{(1)}\right)}{\omega_{\text{MB}}\left(N^{(1)}\right)}\right|$
where $\omega_{\text{MB}}\left(N^{(1)}\right)$ denotes the many-body energy
with the largest distance to
$\omega_{\text{BdG}}\left(N^{(1)}\right)=N^{(1)}\cdot\text{BdG(1)}$, grows
quickly to more than $10\%$. It is emphasized again that the BEC is only
slightly depleted by approximately $1\%$. For the fragmented BEC with
$\Lambda=2.0$ in the deep triple well [panel (d)], the differences become even
more substantial because already for the single-particle excitation
($N^{(1)}=1$) the relative error is about $5\%$. This means that in this case
BdG does not only miss all multi-particle excited states in the low-energy
spectrum, but it is also very inaccurate for the single-particle excitations.

Fig. 24 shows a similar analysis for $N=10$ bosons in a larger shallow lattice
with 10 sites. In principle, the described many-body effects observed in the
triple well also appear in this extended lattice. However, with respect to the
degree of fragmentation, which is $f=0.84\%$ for $\Lambda=1.0$, one can deduce
that these many-body effects set in even earlier than in the triple well
because a comparable depletion is achieved already for weaker repulsion.
Hence, an accurate many-body description is of even larger importance if the
lattice size is increased.

Figure 26: (a) Energies of the zero-quasimomentum modes (ZQMs) for the system
in Fig. 23(a) with $\Lambda=4.0$, up to the top of the second single-particle
band [denoted by BdG(3)]. All shown states are not obtained by the BdG mean-
field approach. The ZQMs are categorized by the symmetry of their density
responses, which is either gerade or ungerade on each lattice site. Mind the
different index $n_{z}$, which enumerates only the ZQMs compared to the index
$n$ of the upper figures which enumerates all excitations. (b) and (c): Real
part of the density responses of the states $(3,0;0,0)$ and $(0,3;0,0)$ [panel
(c)] as well as $(4,1;0,0)$ and $(1,4;0,0)$ [panel (d)]. The vertical lines
separate the lattice sites from each other. Notice the different scales. All
quantities are dimensionless. See text for more details. The figures are taken
from Ref. [168].

The results shown in Fig. 25 demonstrate numerical convergence with respect to
the number of self-consistent orbitals $M$ utilized. In all cases, $M=7$
orbitals clearly suffice for the low-energy excitations. In particular, for
the middle panel describing a BEC in the deep triple well, one can see that
firstly, $M=3$ orbitals already give very accurate values for the excitation
energies, and, secondly, all multi-particle states build up by the states
$(+1)$ and $(-1)$ from the first single-particle band are obtained (65 in
total). To remind the reader, it is stressed again that the BEC is fragmented
in this case, and nevertheless LR-MCTDHB is capable to give converged results.

So far, it was shown that the BdG energies for a BEC in a 1D lattice can be
largely inaccurate, even for weakly-depleted condensates with $f\approx 1\%$.
Though it has been shown analytically that for a BEC with sufficiently weak
and long-ranged repulsion in the Hartree limit, i.e., where
$N\rightarrow\infty$ with fixed $\Lambda$, the BdG equation yields the exact
energies of all single-particle excitations, and the energies of the multi-
particle states are given by multiples and sums of the BdG energies [185]
(extending the findings of a previous work on homogeneous systems [186]).
Whether this result is applicable also for BECs with contact interaction
remains unclear. However, naturally the question arises how fast, if at all,
the excitation spectrum of a BEC with arbitrary shape of the repulsion
approaches the BdG energies upon increasing the number of bosons. This
question is discussed further in the subsequent section, Section 6.3.3, where
excited states in rotating BECs are analyzed.

Figure 27: Low-energy Fourier spectrum of the time evolution of the position
variance per particle in $x$-direction, $\frac{1}{N}\Delta^{2}_{\hat{X}}$, due
to a quantum quench of the lattice depth from $V_{0}=4.9$ to $V_{0}=5.0$. The
system parameters are the same as in Fig. 23(a) with $\Lambda=4.0$. The many-
body spectrum (dashed red), obtained from the evolution with $M=5$ orbitals,
shows two sharp peaks referring to two-particle excitations and one less
intense peak referring to a three-particle excitation. The first intense peak
refers to the two-particle excitation $(1,1;0,0)$, whereas the second intense
peak refers to $(1,0;0,1)$ [or $(0,1;1,0)$]. The less intense peak refers to
$(3,0;0,0)$ [or $(0,3;0,0)$]. All obtained peaks are ZQMs with gerade
symmetry. The GP spectrum (solid black) in the shown energy range is
completely flat. The depletion of the BEC is $f\sim 1.1\%$ during the entire
propagation time. All quantities are dimensionless. See text for more details.

Finally, the zero-quasimomentum modes (ZQMs) are briefly discussed. The latter
are defined in this work as the excitations where $\text{MOD}(P,L)=0$, with
$P$ denoting the total quasimomentum up to the second single-particle band
given by

$P=n_{+1}^{(1)}+n_{+1}^{(2)}-n_{-1}^{(1)}-n_{-1}^{(2)},$ (192)

and $L$ denoting the number of lattice sites. The spectrum of the ZQMs up to
BdG(3), which marks the top of the second single-particle band, is shown in
Fig. 26(a). The system parameters are the same as in Fig. 23(a) with
$\Lambda=4.0$. It is worth noting that all of the shown ZQMs are multi-
particle excitations, meaning that BdG does not account for a single ZQM in
the low-energy spectrum although there exist more than 20. By analyzing the
response densities of the ZQMs, it was found that they can be categorized by
their symmetry with respect to the individual wells, which is either gerade or
ungerade. Specific examples are given in Figs. 26(b) and (c). The ZQMs are
most likely the easiest states to excite in an experiment where the BEC is
subject to either a quench of the potential or the interaction strength. For
example, a simple quench scenario to excite the gerade ZQMs would be a sudden
change of the lattice depth. To this end, the ground state of a BEC with
$N=10$ bosons in the triple well with interaction parameter $\Lambda=4.0$ and
lattice depth $V_{0}=4.9$ is calculated. Then, the depth is suddenly quenched
to $V_{0}=5.0$, and the subsequent time-evolution is investigated. Fig. 27
shows the Fourier spectrum of the dynamical evolution of the position variance
per particle, $\frac{1}{N}\Delta^{2}_{\hat{X}}$, obtained from both
MCTDHB($5$) and GP. Shown is the energy range $5.0\leq\omega\leq 15.5$. At the
many-body level, 3 distinct peaks are observed, whose positions correspond to
the first, second, and fourth gerade excitations of the spectrum in Fig.
26(a). They can thus be labeled according to the occupation of (+1) and (-1)
from the first and second single-particle band. The peaks corresponding to the
two-particle excitations at $\omega\sim 6.35$ and $\omega\sim 13.48$ are more
intense than the peak corresponding to the three-particle excitation at
$\omega\sim 9.02$. Excitations where higher particle numbers are involved are
not intense enough to be distinguished from the background noise. At the GP
level, the spectrum is flat, which means that none of the ZQMs obtained at the
many-body level is included in the mean-field dynamics, which was anticipated
due to the multi-particle nature of the ZQMs.

To summarize, it was demonstrated that the excitation spectrum of a weakly-
depleted BEC in a 1D optical lattice shows clear many-body effects which are
not accounted for at the mean-field level. A depletion of $f\sim 1\%$ is
sufficient to observe these effects. In larger lattices, the sensitivity to
many-body features even increases. It was further shown that a special type of
many-body excited states (ZQMs) are easily accessible by performing a simple
quantum quench to the system. A promising quantity to detect these states is
the position variance, as already seen in the previous section.

####  6.3.3 Rotating Bose-Einstein condensates in an anharmonic trap

After the discussion of many-body excited states in 1D trapped BECs, an
example in 2D, namely the case of rotating BECs in an anharmonic, radially-
symmetric potential crater, is examined. Of particular interest is the weakly-
interacting regime in which the ground-state depletion is marginal and one
would naively expect mean-field approaches like the BdG theory to be valid and
accurate. The results below were recently published in Ref. [187] and it is
referred to this work for additional details.

Rotating BECs were of high interest in the past. Research dealt with the
occurrence of quantized vortices [161, 188, 162, 163, 164] or with the
similarity to the fractional quantum Hall effect for charged particles in the
vicinity of a magnetic field [189, 191, 190, 192]. Moreover, low-lying
excitations in such systems were addressed, for example in vortex lattices
[26, 193, 194, 195]. Furthermore, investigations were made on the decay of the
counter-rotating quadrupole mode [196], on Tkachenko modes [197, 198], on the
twiston spectrum [199], or on excitations in anharmonic traps [200, 201]. In
the latter works, the mean-field excitations were calculated by utilizing the
BdG theory. Computations employing a many-body description arerather rare. An
examples is the analysis of the yrast spectra in a harmonic confinement
obtained by exact diagonalization using the lowest Landau level approximation
[203, 204, 202, 205].

As mentioned above, accurate many-body results for the low-energy spectrum of
rotating BECs in an anharmonic, radially-symmetric trapping potential are
presented. Especially, the effect of the angular momentum in the BECs’ ground
states on the excitation energies is discussed. As for the examples in the
earlier sections, the BdG and many-body results are compared to each other.

Figure 28: (a)-(e) Ground-state densities of rotating BECs with $N=10$ bosons
and interaction parameter $\Lambda=0.5$ for different vorticities $l$. Results
are shown for MCTDHB($7$), the corresponding GP densities look alike (not
shown). For non-zero vorticity, the ground state is a single vortex whose core
size grows with its charge $l$. (f) Low-energy excitation spectra atop the
ground states from the upper panels. Many-body results obtained by LR-
MCTDHB($7$) are denoted by red and blue squares, the BdG results are denoted
by black and blue lines. Once the ground state has non-zero vorticity,
deviations between the mean-field and many-body excitation energies are
observed, and they grow with $l$. Moreover, the density of (many-body) excited
states grows strongly. All ground states are only slightly depleted, with a
maximum of $f\approx 0.38\%$ for $l=4$. Symbols colored in blue refer to the
(+1) excitation energies described in the main text. All quantities are
dimensionless. See text for details. The figure is adapted from Ref. [187].

The crater in which the BEC is confined reads

$V(r)=\begin{cases}C\,e^{-0.5\,(r-R_{C})^{4}}\,\,&\mbox{if }r\leq R_{C}\\\
C\,\,&\mbox{if }r>R_{C}\end{cases}$ (193)

which is essentially the same as in Eq. (179) but without the central ring-
shaped barrier. Here, the radius of the crater is set to $R_{C}=3.0$ and the
height of the wall is kept at $C=200.0$. Furthermore, the same Gaussian
repulsion as in Eq. (180) is employed.

Excitations in this system are analyzed in the rotating frame of reference,
where the Hamiltonian $\hat{H}$ in Eq. (2) contains an additional contribution
accounting for the rotation around the $z$-axis. The whole Hamiltonian thus
reads

$\hat{H}_{\text{rot}}=\hat{H}-\Omega_{\text{rot}}\hat{L}_{z}$ (194)

where $\Omega_{\text{rot}}$ denotes the angular velocity around the $z$-axis
and $\hat{L}_{z}$ denotes the operator of the angular momentum in the
$z$-direction.

A LR analysis is applied atop the ground states of $N=10$ bosons with weak
repulsion strength $\Lambda=0.5$ for different values of
$\Omega_{\text{rot}}$, and thus of the angular momentum per particle $l$,
which is referred to as the vorticity or the charge below. The degree of
condensation, $N(1-f)$, ranges from $9.999$ ($l=0$) to $9.962$ ($l=4$) out of
10 particles, meaning that the BEC is highly condensed for all considered
parameter values. The upper panels in Fig. 28 show the ground-state densities
calculated with MCTDHB($7$) for vorticities $0\leq l\leq 4$. Whereas for
$l=0$, the density resembles a Gaussian in the center of the crater, the shape
of the densities for $l>0$ indicate that the ground state is a vortex with a
density node at its core. One can further see that this core is growing with
$l$. The corresponding GP densities look alike (not shown).

The lower panel in Fig. 28 compares the low-energy excitation spectrum
calculated atop these ground states, both for the mean-field BdG case (black
and blue lines) and for LR-MCTDHB($7$) (red and blue squares). The excitations
denoted by blue symbols are explained below. At first, the spectrum for $l=0$
is examined. Counting the number of excited states, one obtains 2 states from
the BdG approach and 3 states from the many-body approach. Moreover, the
energies of the two mean-field excitations coincide with the corresponding
many-body results. Since the BdG theory only accounts for single-particle
excitations, the third excited stated obtained from LR-MCTDHB is an excitation
where more than one particle is excited from the ground state. In particular,
it represents the state where two bosons are put into the orbital
corresponding to the first excitation. The latter single-particle excitation
is referred to as (+1) since it describes the case of putting a boson from the
ground-state mode with angular momentum $l_{z}=0$ to an orbital with
$l_{z}=1$. Examining the spectra for $l>0$, one observes two major changes. At
first, the density of excitations in the low-energy spectrum clearly grows
with $l$, which can be deduced from the large number of many-body excited
states entering. The number of mean-field excitations also grows (from 2
states for $l=0$ to 4 states for $l=4$), but very slightly compared to the
number of many-body states. Therefore, the amount of excited states that BdG
misses is strongly growing with the vorticity of the BEC’s ground state. The
second observation is that the values of the single-particle energies from BdG
and LR-MCTDHB do not coincide any longer, as it was the case for $l=0$. In
fact, the entire BdG spectrum for vorticity $l=4$ appears to be completely
unrelated to the corresponding many-body spectrum. The increasing gap between
the (+1) results of the two approaches is clearly visible.

Fig. 29 addresses the evolution of the energy of (+1), denoted by
$\omega_{(+1)}$, with respect to the angular velocity $\Omega_{\text{rot}}$.
Shown are the results obtained from the BdG theory, as well as many-body
results for different numbers of bosons ($N=10,\,100$ and $1000$). In the
$l=0$ sector, all energies coincide, indicating that for this case the BdG
theory is sufficient. Henceforth, this obviously changes. With growing vortex
charge of the ground state, the mean-field and many-body lines for (+1)
separate further from each other. Most importantly, increasing the number of
bosons in the condensate only marginally changes this, meaning that the many-
body and mean-field results for large $N$ approach each other very slowly (see
inset). As mentioned in the previous section, it was proven analytically that
for a trapped BEC with sufficiently weak and long-ranged repulsion, the
excitation energies in the Hartree limit converge towards the BdG energies
[185]. The same is believed to happen for a rotating BEC under certain
conditions [206].

Figure 29: Excitation energies $\omega_{(+1)}$ for different angular
velocities $\Omega_{\text{rot}}$. Compared are the BdG energies and the many-
body results for different numbers of bosons $N$, obtained by LR-MCTDHB. The
interaction strength is $\Lambda=0.5$. The transition from ground-state
vorticity $l$ to $l+1$ between two adjacent analyzed values of
$\Omega_{\text{rot}}$ are denoted by dotted vertical lines (mean-field in
black and many-body in red). The higher the ground-state vorticity $l$, the
larger is the separation between the BdG and many-body energies. Inset:
Magnified view for $\Omega_{\text{rot}}=3.8$. All quantities are
dimensionless. See text for details. The figure is adapted from Ref. [187].
Figure 30: Evolution of the energetic difference $\Delta_{N}$ between the
mean-field and many-body results for $\omega_{(+1)}$ for different particle
numbers $N$. The differences are shown relative to $\Delta_{10}$, i.e., the
gap for $N=10$ bosons. Solid lines denote exponential fits. The gap becomes
larger as $l$ grows. Even for $l=1$ and $N=1000$, it is approximately $97\%$
of $\Delta_{10}$, and the slopes of the fit curves do not indicate a sharp
decent towards zero if $N$ is increased further. All quantities are
dimensionless. See text for details. The figure is taken from Ref. [187].

In Fig. 30, the dependence of the difference $\Delta_{N}$ between the mean-
field and many-body energies of (+1) on the number of bosons in the BEC is
analyzed in more detail. The results are given relative to the energy for
$N=10$ bosons, denoted by $\Delta_{10}$. For all shown vorticities, the gap
becomes smaller with increasing $N$. However, it is obvious that even for BECs
with an experimentally relevant number of bosons ($N=1000$) one is still far
away from the predictions of the BdG theory. For $l=1$, one obtains that the
gap size $\Delta_{1000}$ is still about 97$\%$ of $\Delta_{10}$, and
correspondingly higher for $l>1$. Additionally, the shape of the exponential
fit curves to the numerical data does not show any sign of a steep descend
when $N$ is increased further. It can firstly be deduced that the excitation
energies approach the BdG predictions for $l>0$ very slowly, and, secondly,
one apparently needs a very high number of particles in order to enter the
regime where BdG yields approximately the correct results. However, it is
important to note that the convergence towards the BdG excitation energies in
the Hartree limit is not proven for the system under consideration, since Ref.
[206] only deals with the case of a rotating BEC with broken radial symmetry.
The latter situation appears, e.g., for a BEC in a harmonic confinement when
the rotation is fast enough to generate multiple quantized vortices in the
ground state (see, e.g., Ref. [207]). Although the numerical results cannot
give a definite answer for the Hartree limit, they at least suggest that the
regime in which one is clearly in need of an accurate many-body theory for the
excitation energies covers the experimentally relevant regime of rotating
BECs, with $N$ being of the order of $10^{3}-10^{6}$ bosons. A further
discussion on the applicability of the GP and BdG theories in the Hartree
limit can be found in Appendix A.

Figure 31: Low-energy Fourier spectrum of the time evolution of the position
variance in $x$-direction, $\frac{1}{N}\Delta^{2}_{\hat{X}}$, due to a quantum
quench of the radius $R_{C}$ from $3.01$ to $3.0$. The number of bosons is
$N=10$. The many-body spectrum (red), obtained from the evolution with $M=7$
orbitals, shows two sharp peaks in the considered energy range, whereas the GP
spectrum (black) is flat. The first peak refers to the two-particle excitation
$(1_{+1},1_{-1})$, whereas the second peak indicates a de-excitation between
two higher-lying states (see right inset). The left inset shows the evolution
of $\frac{1}{N}\Delta^{2}_{\hat{X}}$ from $0\leq t\leq 30$. The depletion of
the BEC stays below $0.4\%$ during the entire propagation. All quantities are
dimensionless. See text for details.

It is finally suggested, similar to the 1D systems discussed in Section 6.3.1,
how a quench can be utilized in order to excite breathing modes. To this end,
the ground state of a BEC with $N=10$ bosons and angular velocity
$\Omega_{\text{rot}}=3.2$ in the same trap as described in Eq. (193), but with
a slightly larger radius of $R_{C}=3.01$, is calculated. The obtained ground
state has the same symmetry and vorticity $l=3$ as the one for $R_{C}=3.0$ and
$\Omega_{\text{rot}}=3.2$. Then, the crater is suddenly shrinked back to
$R_{C}=3.0$ and the dynamical evolution is computed with MCTDHB($7$), which is
sufficient in order to obtain converged dynamics for the position variance. As
could be seen in the previous sections, the latter quantity is very sensitive
for many-body effects, and therefore, a spectral analysis of its dynamics by
computing the Fourier spectrum of $\frac{1}{N}\Delta^{2}_{\hat{X}}(t)$ is
carried out. Due to the radial symmetry, one could as well analyze the
position variance in the $y$-direction. The result are shown in Fig. 31. The
many-body spectrum shows two distinct peaks, where the first corresponds to
the $(1_{+1},1_{-1})$ two-particle excitation, where one boson is excited to
an orbital with $l=4$ and one is excited to an orbital with $l=2$, such that
the total angular momentum is unchanged. The energy that LR-MCTDHB($7$)
predicts for this state is $\omega=0.605$ which coincides with the peak
position of the Fourier analysis. The second peak at $\omega\approx 1.12$
denotes a de-excitation between two states that are much higher in energy
($\omega\approx 10.6\rightarrow 9.47$, see right inset in Fig. 31). The
occurrence of de-excitations is an indication that the quench described above
already leads to second-order effects, meaning that one has left the LR
regime. With regard to the mean-field spectrum from the GP dynamics, one
observes that in the shown energy range, the Fourier spectrum is completely
flat. The absence of the first peak is anticipated because the GP theory does
not account for two-particle excitations by construction. The absence of the
second peak can be understood from the right inset where the Fourier spectrum
in the interval $9\leq\omega\leq 11$ is shown. One sees that in the GP case,
there is only one peak at $\omega\approx 10.6$, but no second one at
$\omega\approx 9.47$ as in the many-body spectrum. Thus, very different
spectra from the dynamics of the variance at the mean-field and many-body
levels are obtained, and it is shown that the suggested quench leads to an
excitation of a many-body breathing mode. It is stressed that the depletion
stays below $0.4\%$ for the entire propagation, which means that the BEC is
essentially condensed.

To summarize, the low-energy spectrum of a rotating BEC in an anharmonic,
radially-symmetric crater shows significant many-body effects, even for
weakly-interacting, highly-condensed systems. These effects can be observed
already for the lowest-in-energy, single-particle excited states, and remain
even for mesoscopic and large BECs. Furthermore, a simple quench triggers a
many-body breathing mode, which can be observed, e.g., by analyzing the
Fourier spectrum of the position variance.

##  7 Summary and Conclusions

This review presents recently developed many-body approaches capable to
describe and compute accurately the ground state and dynamics as well as the
energy spectrum of trapped ultracold bosons. It is relevant to employ methods
which are on the same footing. The dynamics is calculated by the well-known
MCTDHB method, the ground state is obtained by imaginary-time propagation
using the same method, and the energy spectrum of the system is then
determined by applying the linear-response theory to the ground state computed
via MCTDHB. The resulting theory is named LR-MCTDHB.

The MCTDHB theory is briefly discussed and more attention is paid to the less
spread LR-MCTDHB theory. As the implementation of the MCTDHB is well
documented in the literature, only references are given here, while the newly
developed LR-MCTDHB implementation which makes the theory widely applicable,
is presented here in some detail for the first time. The latter posed a major
challenge due to both the complexity of the underlying theory and the
accompanying technical difficulties like the vast memory consumption and the
necessity of parallelization. The code structure is explained in detail both
with respect to the construction of the linear-response matrix and its
subsequent partial diagonalization which provides the desired excitation
energies. A comprehensive description of the procedure used to calculate the
lowest eigenvalues is also provided. The corresponding benchmarks of the LR-
MCTDHB code (see also Appendix C) with an analytically solvable many-boson
non-trivial model demonstrate its correct functioning and show that the theory
provides numerically exact results.

The second main goal of this review is to present applications to many-body
dynamics and excitations of BECs in interesting scenarios, and, of course, to
discuss the underlying physics. Several systems of trapped BECs in one and two
dimensions are addressed in detail and the main attention is drawn to the
detection of many-body effects that are not describable by applying mean-field
theory. Particularly emphasized are the tunneling dynamics in double-well
systems and the development of fragmentation as a function of time in
initially coherent BECs. The applications shown in this work have demonstrated
that a full many-body treatment of the dynamics of trapped interacting BECs is
inevitable. For instance, it is exhibited that the out-of-equilibrium
tunneling dynamics of initially coherent BECs in 1D and 2D double wells are of
many-body nature. After a few Rabi cycles, tunneling between the wells becomes
suppressed at the many-body level, whereas unperturbed oscillations are
predicted at the mean-field level. Surprisingly, the degree of fragmentation
in 1D is universal for constant values of the interaction parameter. Going to
2D, the sensitivity of the tunneling dynamics to many-body effects is even
enhanced as soon as the bosons carry angular momentum. This can also be
deduced from the larger amount of orbitals necessary in the limit of weak
repulsion as compared to the case of bosons with zero angular momentum.
Further applications deal with the nucleation of phantom vortices, which are
pure many-body objects not observable in the condensate density, as well as
with the tunneling of bosons to open space. In both cases, the initially
coherent systems become fragmented and are, therefore, in need of an accurate
many-body description.

With respect to excitations, it is observed that the ground-state depletion
has a strong impact on the low-energy spectrum. Even for marginal depletion
where one might expect mean-field approaches to give accurate results,
substantial many-body effects like the shift of excitation energies and the
appearance of complicated many-boson excitations occur. Applications in 1D
explores the many-body effects in double- and triple-well systems as well as
in larger lattices. In 2D, it is found that once the bosons carry angular
momentum and the ground state becomes a vortex, the deviations between the
mean-field and many-body excitation spectra grow strongly. Moreover, this
discrepancy tends to persist when the number of particles in the condensate
increases up to an experimentally relevant size. This has a potentially large
impact on experiments on ultracold bosons since the dynamics are mostly
determined by the lowest-in-energy excitations. Hopefully, these findings and
the new insights into the physics of ultracold bosons will stimulate future
experiments.

Before closing, we would like to add several remarks concerning further
general and technical developments of the methods which can further broaden
their range of applicability. Until now the applications of the many-body
linear-response theory were for excitations on top of the system’s ground
state. We emphasize that the LR-MCTDHB can also be applied to any other state
obtained via MCTDHB. In other words, it is possible to compute excitations on
top of any other state. This will help to address excitations which are
difficult to reach from the ground state. Another general possibility lies in