# Probing the edge between integrability and quantum chaos in interacting few-
atom systems

Thomás Fogarty  Quantum Systems Unit, Okinawa Institute of Science and
Technology Graduate University, Onna, Okinawa 904-0495, Japan  Miguel Ángel
García-March  Instituto Universitario de Matemática Pura y Aplicada,
Universitat Politècnica de València, E-46022 València, Spain  Lea F. Santos
Department of Physics, Yeshiva University, New York, New York 10016, USA  N.L.
Harshman  Department of Physics, American University, 4400 Massachusetts Ave.
NW, Washington, DC 20016, USA

###### Abstract

Interacting quantum systems in the chaotic domain are at the core of various
ongoing studies of many-body physics, ranging from the scrambling of quantum
information to the onset of thermalization. We propose a minimum model for
chaos that can be experimentally realized with cold atoms trapped in one-
dimensional multi-well potentials. We explore the emergence of chaos as the
number of particles is increased, starting with as few as two, and as the
number of wells is increased, ranging from a double well to a multi-well
Kronig-Penney-like system. In this way, we illuminate the narrow boundary
between integrability and chaos in a highly tunable few-body system. We show
that the competition between the particle interactions and the periodic
structure of the confining potential reveals subtle indications of quantum
chaos for 3 particles, while for 4 particles stronger signatures are seen. The
analysis is performed for bosonic particles and could also be extended to
distinguishable fermions.

##  1 Introduction

The interest in quantum chaos, especially when caused by the interactions
between particles, has grown significantly in the last few years due to its
relationship with several questions of current experimental and theoretical
research that arise in atomic, molecular, optical, condensed matter, and high
energy physics, as well as in quantum information science. In interacting
many-body quantum systems, quantum chaos ensures thermalization and the
validity of the eigenstate thermalization hypothesis (ETH) [114, 17, 28],
hinders localization [98, 79, 74], facilitates the phenomenon of many-body
quantum scarring [8, 107], leads to diffusive transport [12, 55] and causes
the fast spread of quantum information [91, 75].

Quantum chaos refers to properties of the spectrum and eigenstates that appear
in the quantum domain when the classical counterpart of the system is chaotic
in the sense of mixing and positive Lyapunov exponent. The features are
similar to what we find in random matrix theory [76], namely the eigenvalues
are strongly correlated [44] and the eigenstates in the mean-field basis are
close to random vectors [114]. This quantum-classical correspondence is well
established for systems with few degrees of freedom [103], such as billiards,
the kicked rotor, and the Dicke model, where the source of chaos is
respectively the shape of the billiard, the strength of the kicks, and the
collective interaction between light and matter. In the case of systems with
many interacting particles, the semiclassical analysis is challenging and
sometimes not well defined, so the common approach has been to refer to many-
body quantum systems that present the above mentioned properties of spectrum
and eigenstates as chaotic, even when the classical limit is not analyzed.

The purpose of this work is to identify a minimum model of interacting
particles that is chaotic and that can be experimentally studied with cold
atoms. There are theoretical examples in the literature of quantum systems
with only 3 or 4 interacting particles that already exhibit chaotic
properties. They include the cesium atom, which has 4 valence electrons [33];
systems composed of 4 particles of unequal masses in a harmonic trap [49] and
3 particles with unequal masses on a ring [54]; 3 or 4 excitations in spin-1/2
chains with short-range [94] or long-range couplings [115]; and even spin-1/2
chains with only 3 sites [77]. In the context of thermalization due to chaos,
we also find works that obtained the Fermi-Dirac distribution in systems with
only 4 particles [97, 31, 34, 32, 57].

We consider a one-dimensional system with $N$ identical particles that is
split into wells separated by delta-function barriers. A single barrier
defines a double-well system and many barriers results in the finite
Kronig–Penney model [65, 89]. We focus on the sector of states symmetric under
particle exchange and spatial parity, and that the particles interact via
contact interactions which are modeled with a delta function as in the Lieb-
Liniger Hamiltonian [72, 81]. The system is integrable when there are finite
barriers _or_ finite interactions (not both). It is also solvable in the
limiting cases of infinite barrier strength and infinite interaction strength.
However, we provide numerical evidence that when the interaction strength and
the barrier strength are simultaneously finite, integrability is broken. In
this case, strong signatures of quantum chaos emerge for $N=4$ particles in
the presence of just one barrier (double-well system). We also demonstrate
that the signatures of quantum chaos get enhanced as we increase the number of
barriers, in which case strong level repulsion is verified for as few as $N=3$
particles. Our analysis is done for bosons, but can also be extended to
systems with a small number of distinguishable fermions.

##  2 Experimental realization

The experimental realization of our model can be done with a controllable
number of interacting atoms trapped in one or several one-dimensional traps.
Like experiments based on atoms in optical lattices [14, 70], such systems
allow precision control and find potential applications in quantum engineering
and quantum technologies [102]. They are also suitable test beds for
addressing the question of the transition from few- to many-body systems [16,
15, 102]. Our choice of the Kronig-Penney potential is also motivated by
recent experiments that have achieved coherent optical lattices with a sub-
wavelength spatial structure that consists of ultra-narrow barriers [109].
This technique can create sharp-box potentials and the heights of individual
barriers can be varied, which adds a further tunable parameter to the system.

Our results apply to few-body atomic systems with identical particles that
possess spatial wave functions symmetric under particle exchange. Wave
functions with this spatial symmetry are clearly relevant for the description
of bosons, but our results are equally valid for distinguishable fermions,
i.e. identical fermions with internal degrees of freedom like spin. For a
system of distinguishable fermions, the antisymmetry required by particle
statistics can be carried by the spin or internal wave function. All permitted
symmetries for three fermions or bosons (distinguishable or indistinguishable)
are discussed in [46, 36, 50]; for more particles in [47]; for few particles
in double or few wells in [37, 48].

Due to losses via three-body re-combination, experiments with a few bosons
trapped in the ultracold regime have an additional difficulty when compared
with those with a few fermions. However, for the order of tens of bosons, it
was shown in [18] that one can successfully load dipole traps by means of
evaporative cooling. Smaller number of atoms can be loaded in optical lattices
[112], in arrays of double wells [25], or in a two-site optical ring, which
can be appropriately reshaped into a Gaussian trap [52]. An alternative
experimental route leading to cooled atoms trapped in several wells is that of
few atoms in optical tweezers, which for a single atom in the ground state was
accomplished in [60]. In subsequent papers it was experimentally demonstrated
the trapping of two ${}^{87}\mathrm{Rb}$ bosonic atoms in two wells [61] or in
uniformly filled arrays of traps [68]. Recent advances allow the laser cooling
of atoms in optical tweezers [88]; the trapping of individual atoms in optical
tweezer arrays [93]; and even the loading of atoms one by one [99, 4, 30] in a
one-dimensional array [30].

Another experimental context for the results in this paper is the ground
breaking experiments that showed the accessibility and versatility of systems
with a very small number of interacting fermions. Reference [100] demonstrated
that a deterministic number of ultracold fermions could be extracted from a
larger ensemble by applying a tightly confined one-dimensional dimple
potential. In this experiment, a few ${}^{6}\mathrm{Li}$ atoms in the two
lowest-energy Zeeman substates were trapped. The strength of the interactions
between atoms in different spin states could be controlled via Feshbach
resonance. In a subsequent experiment by the same group, they considered one
atom of one species (an impurity) interacting with an increasing number of
identical fermions, being able to build a small Fermi sea adding fermions one
by one in a controllable manner [110]. One can have more than two components
in a few fermion system, as in the experiment reported in [84], where a one
dimensional system with a tunable number of spin components was realized.

In addition to the controllable number of atoms, another ingredient required
for the realization of our model is the ability to change the trapping
potential, creating double, triple or generally multiwell potentials. The same
group that realized few trapped fermions in [100] was also able to trap few
fermions in one dimensional double-well [78, 7] and multi-well systems [6].

##  3 Model

In the simplest case of equally-spaced barriers, the Hamiltonian describing
our system takes the form

$\frac{1}{\epsilon_{1}}H(N,W,\tau,\gamma)=T^{N}+\tau V^{N,W}+\gamma U^{N}$
(1a) where $\displaystyle T^{N}$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{N}\frac{\partial^{2}}{\partial x_{i}^{2}}$ (1b)
$\displaystyle V^{N,W}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{k=1}^{W-1}\delta(x_{i}-\pi k/W)$ (1c)
$\displaystyle U^{N}$ $\displaystyle=$ $\displaystyle\sum_{\langle
ij\rangle}\delta(x_{i}-x_{j}).$ (1d)

This model realizes a system with $N$ identical interacting particles of mass
$m$ trapped in a one-dimensional box of length $L$ that is disrupted by $W-1$
delta-barriers. In the equation above, $x_{i}\in[0,\pi]$ are the positions of
the particles scaled by the length $L/\pi$, the energy scale is provided by
$\epsilon_{1}=\hbar^{2}\pi^{2}/(2mL^{2})$ (henceforth set to unity), and
$\tau\geq 0$ and $\gamma\geq 0$ are unitless parameters describing the barrier
strength, and interaction strength, respectively. Here we consider repulsive
interactions, as it is the most common scenario for ultracold bosons.

The goal of our analysis of the Hamiltonian $H(N,W,\tau,\gamma)$ is to
understand how the signatures of quantum chaos scale with the number of
particles $N$, the number of wells $W$, the barrier strength $\tau$ and the
interaction strength $\gamma$. The particularly simple form (1) of the
Hamiltonian means that it is amenable to analytic and numeric calculations.

###  3.1 Solvable and integrable limiting cases

For a fixed $N$ and $W$, consider the parameter space window
$\tau\in[0,\infty)$ and $\gamma\in[0,\infty)$, depicted in Fig. 1. An
interesting feature of the Hamiltonian $H(N,W,\tau,\gamma)$ is that the model
has exact solutions for four special limiting cases:

  1. 1.

$H(N,W,0,0)=T^{N}$: $N$ non-interacting identical particles in a one-
dimensional infinite square well with width $L$.

  2. 2.

$H(N,W,0,\infty)$: $N$ hard-core identical particles in a one-dimensional
infinite square well with width $L$.

  3. 3.

$H(N,W,\infty,0)$: $N$ non-interacting identical particles distributed in $W$
identical one-dimensional infinite square wells with width $L/W$.

  4. 4.

$H(N,W,\infty,\infty)$: $N$ hard-core identical particles distributed in $W$
identical one-dimensional infinite square wells with width $L/W$.

In all of these four cases, the configuration space is sectioned into one or
more $N$-dimensional polytopes with high symmetry. Solving for the spectrum is
equivalent to solving the Schrödinger equation with Dirichlet boundary
conditions. Exact solutions for the Schrödinger equation in these polytopes
can be constructed from symmetrized combinations of one-particle states using
methods of Refs. [64, 108, 59, 82].

Figure 1: Depiction of the (compactified) $(\tau,\gamma)$ parameter space
window $\tau\in[0,\infty)$ and $\gamma\in[0,\infty)$. At the four corners of
the window, the Hamiltonian $H(N,W,\tau,\gamma)$ has exact solutions and is
superintegrable. The four edges connecting these corners are integrable models
where exact solutions can be found by the solution of coupled transcendental
equations. Except for possibly the case of $N=2$, the $H(N,W,\tau,\gamma)$
does not appear integrable or solvable for arbitrary $(\tau,\gamma)$ in this
parameter space

Beyond these four exactly solvable special cases, the Hamiltonian
$H(N,W,\tau,\gamma)$ is also _integrable_ along the four edges of the
$(\tau,\gamma)$ parameter space window (cf. Fig. 1). Referring to the four
‘corner’ models denoted above, the integrable limits are:

  * $1\leftrightarrow 2$

$H(N,W,0,\gamma)$: Without barriers, the Hamiltonian (1) for $N$ identical
particles in an infinite square well with delta-function interactions is
solvable by coordinate Bethe ansatz [39, 40, 80]. The energies are the
solution of coupled transcendental equations that depend on $N$ and $\gamma$.

  * $1\leftrightarrow 3$

$H(N,W,\tau,0)$: $N$ non-interacting particles in an infinite square well with
$W-1$ barriers. The Hamiltonian separates into $N$ identical one-dimensional
sub-Hamiltonians. The spectrum for each sub-Hamiltonian is obtained via
solution to transcendental equations that depend on $\tau$ [89].

  * $2\leftrightarrow 4$

$H(N,W,\tau,\infty)$: In this limit of finite wells with hard-core
interactions, the solutions are Tonks-Girardeau constructions derived from the
symmetrized Slater determinants of the solutions of the non-interacting case
$H(N,W,\tau,0)$ [41].

  * $3\leftrightarrow 4$

$H(N,W,\infty,\gamma)$: With infinite barriers, the configuration space
fractures into $W^{N}$ configurations, i.e. each particle
$i\in\\{1,\ldots,N\\}$ is in a well $j\in\\{1,\ldots,W\\}$ with width $L/W$.
Each of these particle-well configurations is an independent system, since
with $\tau\to\infty$ there is no tunneling among wells. Within each
configuration, the solutions are similar to the coordinate Bethe ansatz
solutions of $H(N,W,0,\gamma)$, except that now there may be different
partitions of interacting particles and the size of the well is $L/W$.

Note that these four integrable models allow us to establish an adiabatic map
between the energy levels for the solvable (and superintegrable) models at the
corners of parameter space.

###  3.2 Symmetries and degeneracies

To analyze signatures of chaos, the Hilbert space $\mathcal{H}$ must first be
decomposed into subspaces with fixed symmetry. The Hamiltonian (1) has two
symmetries for any interaction strength and barrier strength in the
$(\tau,\gamma)$ parameter window. First, it is symmetric under particle
permutations $S_{N}$, i.e. any permutation $p\in S_{N}$ is represented as a
linear transformation in configuration space
$(x_{1},x_{2},\ldots,x_{N})\to(x_{p_{1}},x_{p_{2}},\ldots,x_{p_{N}})$ that
leaves the Hamiltonian (1) invariant. The second symmetry is total parity
inversion $\Pi$ about the center of the well, implemented as the affine linear
transformation $x_{i}\to\pi-x_{i}$ for all $i$.

These two symmetries allow eigenstates to be classified by irreducible
representations of $S_{N}$ and of parity $\Pi$ and reduce the total Hilbert
space $\mathcal{H}$ into subspaces with a given symmetry (cf. [47]). In this
work, we focus on the sector of Hilbert space
$\mathcal{H}^{[N]+}\subset\mathcal{H}$ containing states with bosonic symmetry
under particle exchange and positive parity. Note that a special property of
delta-function operators $V^{N,W}$ and $U^{N}$ is that they vanish on certain
subspaces of $\mathcal{H}$, because the support of the delta-functions
coincides with nodal lines of symmetrized eigenstates of $T^{N}$. In
particular, there is a subspace of $\mathcal{H}^{[N]+}$ upon which $V^{N,W}$
vanishes when $N$ is even or both $N$ and $W$ are odd. Also, in the limits
$\tau\to\infty$ and $\gamma\to\infty$, both $V^{N,W}$ and $U^{N}$ must vanish
on any states with finite energy; in this limit the wave functions must have
nodal surfaces that coincide with the support of these operators.

Degeneracies in the spectrum either originate with the symmetries of the
Hamiltonian or they are designated ‘accidental’. Because the only symmetries
for generic $(\tau,\gamma)$ are particle exchanges $S_{N}$ and parity $\Pi$,
there should be no degeneracies originating in symmetry in the interior of
parameter space. This is because the irreducible representations $[N]+$ for
totally symmetric bosonic states are one-dimensional. However, additional
symmetries arise in the limiting cases of zero and infinite strengths $\tau$
or $\gamma$, including the symmetries of separability and system decoupling
due to infinite barriers [48]. These symmetries explain the integrability of
the edge models and the superintegrability of the corner models. When $\tau$
and $\gamma$ are both finite, then these additional symmetries are broken and
we expect these energy levels that cross at certain parameters along the edges
of the model space to repel. Level repulsion, which is a signature of chaos,
is thus expected in the interior of this model space (see Fig. 1), which is
indeed what we verify for $N\geq 3$.

Besides degeneracies originating in symmetries, several other kinds of
accidental degeneracies should be possible in our system. The variation of two
parameters is sufficient to allow for ‘diabolic points,’ topologically-stable
conic degeneracies between two energy levels [9, 11]. These degeneracies can
be distinguished from ‘near misses’ by looking at how the wave function
transforms when small loops in control space are taken and in subsequent work
we will consider such loops. There are also degeneracies for the solvable
corner models arising from the number theory of decomposing integers into sums
of squares of integers. Also called Pythagorean degeneracies [101, 48], the
density of such degeneracies grows slowly with energy but eventually comes to
dominate the spectrum [9]. Finally, there are other degeneracies
characteristic of Bethe-ansatz solvable systems [53] that should be relevant
for some of the integrable edge models. It is not clear how these additional
‘accidental’ degeneracies at the corners and edges affect the spectrum of the
interior of parameter space, but they inform our interpretation of the level
statistics presented below for the $N=2$ case.

###  3.3 Density of states

The density of states of the Hamiltonian $H(N,W,\tau,\gamma)$ is used to
interpret some of the numerical results below. Perhaps surprisingly, it is
independent of $\tau$ and $\gamma$ to leading order in the energy $E$. The
density of states in the sector $\mathcal{H}^{[N]+}$ is

$\displaystyle\rho^{[N]+}(E)$ $\displaystyle=$
$\displaystyle\frac{1}{2^{N+2}(N+1)!}\frac{\pi^{N/2}}{\Gamma(N/2+1)}E^{N/2-1}$
(2) $\displaystyle+O\left[N,W,\tau,\gamma\right](E^{N/2-3/2}),$

where we have used the notation $O\left[N,W,\tau,\gamma\right]$ to indicate
that generally, the coefficient in front of the subleading term proportional
to $E^{N/2-3/2}$ will depend on all the parameters of the Hamiltonian.

Figure 2: Density of states for $N=2$ (blue points), $N=3$ (red points) and
$N=4$ (yellow points). Solid lines show the leading term in Eq. (2) for the
respective particle numbers. Symbols denote different system parameters:
$\gamma=10$ and $\tau=0$ (circles), $\gamma=\tau=10$ and $W=2$ (triangles) and
$\gamma=\tau=10$ and $W=10$ (stars).

We derive Eq. (2) in the Appendix for the four solvable ‘corner’ models, and
infer that it holds for the entire parameter window. More generally, adding
one-dimensional delta-functions to an otherwise free problem should not change
the density of states. By Weyl’s Law [56], the volume of phase space can be
related to the density of states, and boundary conditions (whether Dirichlet,
Neumann or hybrid) do not change the phase space volume to leading order.

In Fig. 2 we show that the leading term in Eq. (2) agrees with the scaling of
the density of states regardless of number of wells $W$, barrier strength
$\tau$ or interaction strength $\gamma$. Because the leading term in Eq. (2)
is the same for the entire parameter window for a given $N$, at least
heuristically we expect that the variation in level statistics with the
parameters $(\tau,\gamma)$ is governed by the properties of the subleading
term (or even lower order terms). In the Appendix, we support this hypothesis
by looking at how energy levels change as the model parameters are
adiabatically tuned along integrable edge models to connect the spectra of the
solvable corner models.

Note that for $N=2$, the density of states is a constant and the subleading
correction is actually decreasing with energy as $E^{-1/2}$. In Fig. 2, all
data points for $N=2$ lie on top of the line given by the leading term in Eq.
(2). For $N=3$, the density of states increases as $E^{1/2}$ and the
subleading term is constant. Only for $N=4$ are both the leading and the
subleading term of the density of states growing with energy. If the
subleading term is important for understanding the density of level crossings
as we propose, then this helps to explain why $N=4$ is the threshold when the
level statistics conform to the expectations of random matrix theory across
such a wide range of parameter space, whereas the $N=3$ only in a limited
range.

##  4 Indicators of chaos

In the parameter window of interest, finite interactions and finite barrier
heights, the Hamiltonian (1) is diagonalized numerically using a basis
consisting of $N$-particle non-interacting eigenstates of $T^{N}$ in the
Hilbert space $\mathcal{H}^{[N]^{+}}$ [46]. In the following we discuss the
different indicators of chaos that are used in this work along with the
numerical results.

###  4.1 Energy spacing distribution

To study the degree of short-range correlations between the eigenvalues, we
use the distribution $P(s)$ of the spacings between neighboring levels
obtained after unfolding the spectrum. In generic integrable models, the level
spacing distribution is Poissonian, $P_{\rm P}(s)=e^{-s}$, when the energy
levels are uncorrelated and not prohibited from crossing, although different
shapes emerge for “picket-fence”-kind of spectra [10, 85, 38] and systems with
an excessive number of degeneracies [113]. For chaotic systems with real and
symmetric Hamiltonian matrices, as the one considered here, the level spacing
distribution follows the Wigner-Dyson distribution [76, 44], $P_{\rm
WD}(s)=(\pi s/2)\exp\left(-\pi s^{2}/4\right)$, which indicates that the
eigenvalues are correlated and repel each other.

First we investigate the minimal case for chaos in our system, where only one
barrier is inserted centrally in the square well, thereby creating a double
well potential ($W=2$). In Fig. 3 we show the level spacing distributions for
$2$, $3$ and $4$ particles for different interactions and at a fixed barrier
height of $\tau=10$.

Figure 3: Energy level spacing distributions for (a,d,g) $N=2$ particles,
(b,e,h) $N=3$ particles and (c,f,i) $N=4$ particles with different interaction
strengths $\gamma$. The barrier height is fixed to $\tau=10$ and the number of
wells is $W=2$. The red solid line is the Poissonian distribution and the blue
solid line is the Wigner-Dyson distribution.

For weak interactions, $\gamma=1$, a picket fence pattern is noticeable for
$N=2$ particles [Fig. 3 (a)], which indicates non-generic correlations in the
energy spectrum [10]. Additionally, the peak at $s=0$ signals excessive
degeneracies typical for solvable systems [113, 53]. For $N=3$, the picket
fence structure vanishes and the distribution is closer to Poissonian with
some remaining evidence of additional degeneracies at $s=0$ [Fig. 3 (b)]. For
$N=4$, the distribution is also close to Poissonian, except for a slight
decrease at $s=0$ that provides evidence for level repulsion already at weak
interactions [Fig. 3 (c)].

At stronger interactions ($\gamma\geq 10$), the distributions for $N=2$
particles in Fig. 3 (d) and Fig. 3 (g) become close to Poissonian and level
repulsion is not evident. This Poissonian distribution suggests that for the
$N=2$ and $W=2$ the Hamiltonian is integrable (or effectively so) for all
values of $(\tau,\gamma)$; see the conclusion for a further discussion on this
point. In contrast, some degree of level repulsion is already noticeable for
$N=3$ [Fig. 3 (e) and Fig. 3 (h)] and the Wigner-Dyson distribution is visible
for $N=4$ particles [Fig. 3 (f) and Fig. 3 (i)].

To explore the crossover of the energy spectrum from Poissonian to Wigner-
Dyson and the role of interactions and the potential barrier, we fit our
numerical results to the Brody distribution [22] (see an alternative in [58]),

$\displaystyle P_{\beta}(s)$ $\displaystyle=$
$\displaystyle(\beta+1)bs^{\beta}\exp(-bs^{\beta+1}),$ (3) $\displaystyle b$
$\displaystyle=$
$\displaystyle\left[\Gamma\left(\frac{\beta+2}{\beta+1}\right)\right]^{\beta+1}.$

For the Wigner-Dyson distribution, $\beta\sim 1$, while the Poissonian
distribution leads to $\beta\sim 0$.

Close to the integrability limits, vanishing interactions ($\gamma\approx 0$)
with finite barriers ($\tau\neq 0$) or finite interactions ($\gamma\neq 0$)
with vanishing barrier height ($\tau\approx 0$), the energy spacing
distributions can display a picket-fence or quasi-Poissonian distribution as
discussed in Fig. 3 for $N=2$, which results in $\beta<0$. We highlight these
non-generic regions as black in Figs. 4 (a)-(f), showing that they are more
prevalent at smaller particle numbers, while their footprint in the parameter
space almost vanishes for $N=4$. As we increase the interactions ($\gamma\gg
0$) and the barrier heights ($\tau\gg 0$), the degeneracies of the $N=2$
system are destroyed and the level spacing distributions become closer to
Poissonian for a large region of the parameter space with $\beta$ remaining
below $0.35$.

Figure 4: Brody distribution parameter $\beta$ as a function of $\tau$ and
$\gamma$ for (a,d) $N=2$, (b,e) $N=3$, and (c,f) $N=4$ particles. The number
of wells are (a-c) $W=2$ and (d-f) $W=10$. For $W=2$ wells the maximum
$\beta_{max}=0.8269$ is found for $N=4$ particles at $\gamma=20$ and
$\tau=12.5$, while for $W=10$, we find $\beta_{max}=0.9512$ for $N=4$
particles at $\gamma=45$ and $\tau=12.5$. Black areas show $\beta<0$
indicating degeneracies and picket-fence spectra.

For more particles, the distributions are different, with indications of
energy level repulsion emerging for $N=3$ in the region
$15\lesssim\gamma\lesssim 45$ and $5\lesssim\tau\lesssim 35$ where $\beta\sim
0.5$ [Fig. 4 (b)]. For $N=4$ the fitting parameter $\beta$ in almost the
entire parameter space is larger than $0.5$, with the areas of integrability
confined to the edges of the $(\tau,\gamma)$ parameter space [Fig. 4 (c)]. In
fact, we find the maximum to be around $\beta\sim 0.83$ at $\gamma=20$ and
$\tau=12.5$, indicating that the energy spacing distribution gets close to
Wigner-Dyson.

In Figs. 4 (d)-(f) we also consider multiple wells, with $W=10$. For more
wells, the density of states does not change to leading order in $\tau$, but a
greater proportion of the eigenstates of $T^{N}$ feel the effect of the
barriers more acutely. For example, for $N=2$ and $W=2$, the barrier potential
$V^{N,W}$ vanishes on half of the eigenstates of $T^{N}$ in
$\mathcal{H}^{[N]+}$, whereas for $N=2$ and $W=10$, the barrier potential
$V^{N,W}$ only vanishes on one tenth of the eigenstates of $T^{N}$ in
$\mathcal{H}^{[N]+}$.

Comparing the case of $W=10$ to $W=2$ in Fig. 4, we see that larger values of
$\beta$ are indeed found for $W=10$ and the region of the parameter space
where $\beta$ is large has increased significantly for both $N=3$ and $N=4$.
This suggests that in systems with a large number of wells, the energy level
repulsion is enhanced and chaos could be observed for systems with as few as
$N=3$ particles. However, for $N=2$, we find a maximum of $\beta\approx 0.5$
and this value does not approach the chaotic limit by increasing the number of
wells. As before, the case of only two particles resists the transition to
chaos. In Sect. 4.4 below, we discuss in more detail how the signatures of
chaos change as the number of barriers is increased.

Figure 5: Inverse of the kurtosis $\mathcal{K}^{-1}_{T}$ as a function of
$\tau$ and $\gamma$ for (a,d) $N=2$ (b,e) $N=3$ and (c,f) $N=4$ particles. The
number of wells are (a-c) $W=2$ and (d-f) $W=10$. For $W=2$ the minimum
kurtosis is $4.3592$ at $\gamma=30$ and $\tau=15$, while for $W=10$ the
minimum kurtosis is $3.3274$ at $\gamma=20$ and $\tau=7.5$. Figure 6: Off-
diagonal elements of the kinetic energy operator $T_{mn}$ for (a-d) $W=2$
wells and (e-h) $W=10$ wells. The probability distributions in the respective
potentials are shown for (a,e) $N=2$ (b,f) $N=3$ and (c,g) $N=4$. (d,h) The
ratio $\Gamma_{T}$ as a function of $\omega=|E_{m}-E_{n}|$ in an energy window
at $E_{\rm{mid}}=700$ and width $\Delta E=100$. The interaction strengths and
barrier heights are chosen at the point of minimum kurtosis for $4$ particles
in Fig.5(c,f), namely $\gamma=30$ and $\tau=15$ for $W=2$, and $\gamma=20$ and
$\tau=7.5$ for $W=10$.

###  4.2 Off-diagonal ETH

Quantum chaos ensures the validity of the ETH, so we can also use the
indicators of ETH to detect the transition to chaos. Two conditions need to be
satisfied for a few-body observable $O$, evolving according to

$\displaystyle O(t)$ $\displaystyle=$ $\displaystyle\sum_{m\neq
n}C_{m}^{*}C_{n}e^{-i(E_{n}-E_{m})t}O_{mn}$ (4)
$\displaystyle{}+\sum_{m}|C_{m}|^{2}O_{mm},$

to reach thermal equilibrium. In the equation above, $C_{m}=\langle
m|\Psi(0)\rangle$ is the overlap between the eigenstate $|m\rangle$ of the
Hamiltonian that describes the system and the initial state $|\Psi(0)\rangle$,
and $O_{mn}=\langle m|\hat{O}|n\rangle$. The first condition is that of
equilibration, which depends on the first term on the right-hand side of Eq.
(4). At long times, due to the lack of degeneracies of chaotic Hamiltonians,
the small values of the coefficients $C_{m}$’s obtained when the systems is
quenched far from equilibrium, and the small values of the off-diagonal
elements $O_{mn}$’s caused by the chaotic eigenstates, the first term in Eq.
(4) leads to small fluctuations that decrease with system size and cancel out
on average. So apart from small fluctuations, the observable reaches its
infinite-time average $\sum_{m}|C_{m}|^{2}O_{mm}$. The second condition is
that this infinite-time average approaches the thermodynamic average as the
system size increases, confirming that the equilibrium is indeed thermal.
These two steps are usually referred to as off-diagonal- and diagonal-ETH,
respectively.

Here, we consider the off-diagonal-ETH to detect the transition to chaos. The
distribution of $O_{mn}$ in chaotic (thermalizing) systems is Gaussian [13,
66, 19, 20, 92, 21], reflecting the chaotic structure of the eigenstates,
while other forms emerge for integrable models [66]. To quantify how close the
distribution of $O_{mn}$ is to a Gaussian, we use the kurtosis,

$\mathcal{K}_{\hat{O}}=\frac{\langle(O_{mn}-\langle
O_{mn}\rangle)^{4}\rangle}{\sigma^{4}},$ (5)

where $\langle.\rangle$ indicates the average over all pairs of eigenstates
$|m\rangle\neq|n\rangle$ and $\sigma$ is the standard deviation of the
distribution. For Gaussian distributions the kurtosis is given by
$\mathcal{K}_{\hat{O}}=3$.

In our calculations of $\mathcal{K}_{\hat{O}}$, we consider the kinetic energy
operator $T^{N}$ (in the following we drop the superscript to simplify the
notation). To ensure no effects from the bottom edge of the spectrum, we
choose a energy window far from the ground state and include only states
within $E_{m}\in[E_{\rm{mid}}-\Delta E,E_{\rm{mid}}+\Delta E]$, where the
center of the energy window is high in the spectrum at $E_{\rm{mid}}=700$ and
its width is $\Delta E=100$. In Figs. 5 (a)-(f) we show the inverse of the
kurtosis, $1/\mathcal{K}_{T}$, with the maximal value $1/3$ indicating a
Gaussian distribution and therefore the presence of chaos.

For $N=4$ particles and $W=2$ [Fig. 5 (c)], the distribution is close to
Gaussian for barrier heights $1\lesssim\tau\lesssim 30$ and interactions
$20\lesssim\gamma\lesssim 60$, which coincides with the region of $\beta>0.8$
in Fig. 4 (c). We choose the minimum value of the kurtosis in Fig. 5 (c) and
show the Gaussian probability distribution of $T_{mn}$ in Fig. 6 (c). In
contrast, for lower particle numbers, $N=2,3$, the distribution of the off-
diagonal elements $T_{mn}$ is sharply peaked, as seen in Fig. 6 (a) and Fig. 6
(b), which is consistent with the kurtosis having values $\mathcal{K}_{T}\gg
3$ in Fig. 5 (a) and Fig. 5 (b).

The kurtosis shows a similar enhancement due to the presence of more barriers,
attaining a minimum of $\mathcal{K}_{T}\approx 3.3$ for $N=4$ particles in
$W=10$ wells [Fig. 5 (f)]. Indeed, for $1\lesssim\tau<20$ and a broad range of
interactions, the kurtosis is close to 3. Increasing the number of barriers
also moves the band of minimal kurtosis to lower values of $\tau$, as lower
barrier heights are necessary to retain the competition with the inter-
particle interactions. In a similar region of the parameter space, there is
also a visible minimum of the kurtosis for $N=3$ particles [Fig. 5 (e)]. In
fact, when taking the interaction and barrier height which give the minimum
value of kurtosis [$\gamma=10$ and $\tau=7.5$ for $N=4$], there are distinct
Gaussian probability distributions for both $N=4$ and $N=3$ particles [see
Figs. 6 (f,g)]. However, for $N=2$ the off-diagonal elements of the kinetic
energy operator $T_{mn}$ do not indicate the presence of chaos in either the
kurtosis [Fig. 5 (d)] or the probability distribution [Fig. 6 (e)].

To quantify how the off-diagonal elements $T_{mn}$ behave as a function of the
energy difference $\omega=|E_{m}-E_{n}|$ we also show in Fig. 6 the ratio

$\Gamma_{T}=\frac{\overline{|T_{mn}|^{2}}}{\overline{|T_{mn}|}^{2}},$ (6)

which is equal to $\pi/2$ for a Gaussian distribution. In Fig. 6 (d) and Fig.
6 (h) this ratio is shown for $W=2$ and $W=10$ wells, respectively, for the
same Hamiltonian parameters that gave the lowest values of the kurtosis
discussed in the preceding paragraph. For $W=2$ wells, the values of
$\Gamma_{T}$ for $N=4$ sit very close to $\pi/2$ over a large range of
$\omega$, confirming that the Gaussianity of the distribution is preserved at
different energy spacings. The ratio for $N=2,3$ has a large variance with the
majority of the points being far from $\pi/2$, which is indicative of the
peaked distributions in Fig. 6 (a) and Fig. 6 (b). For $W=10$ wells, we find
that not only $\Gamma_{T}\approx\pi/2$ for $N=4$ particles, but also the $N=3$
system approaches this result. This suggests that the $N=3$ system can be
tuned between integrability and chaos by changing the trapping potential.

###  4.3 Survival Probability

Spectral correlations get manifested also in the evolution of the survival
probability,

$\displaystyle\langle S_{P}(t)\rangle$ $\displaystyle=$
$\displaystyle\langle\left|\left<\Psi(0)|\Psi(t)\right>\right|^{2}\rangle$ (7)
$\displaystyle=$
$\displaystyle\langle\sum_{n,m}|C_{n}|^{2}|C_{m}|^{2}e^{-i(E_{n}-E_{m})t}\rangle,$

in the form of what is known as the correlation hole [69, 86, 43, 111, 51, 73,
2, 71, 104, 105, 106, 95, 67, 29], recently referred to also as the “ramp”
[27]. The correlation hole corresponds to a dip below $\langle
S_{\infty}\rangle=\langle\sum_{n}|C_{n}|^{4}\rangle$, which is the infinite-
time average (saturation value) of $\langle S_{P}(t)\rangle$. This dip emerges
also in the spectral form factor
$\langle\sum_{n,m}e^{-i(E_{n}-E_{m})t}\rangle$, but contrary to this one, the
survival probability is a true dynamical quantity. In the equation above,
$\left<.\right>$ indicates averages. The survival probability is non-self-
averaging [87, 96], so the correlation hole is not visible unless averages are
performed. They can be done over initial states, disorder realizations, or, as
in our case, they correspond to moving time averages. The correlation hole
detects short- and long-range correlations in the spectrum, and in addition,
it does not require unfolding the spectrum or separating it by symmetries [29,
92]. In cold atom systems the survival probability is commonly used to probe
the non-equilibrium dynamics of few- [23, 62] and many-body systems [42, 63,
26, 35], and can be experimentally measured using interferometric techniques
[24].

Figure 7: Survival probability of (a,d) $2$ particles (b,e), $3$ particles,
and (c,f) $4$ particles. The number of wells is (a-c) $W=2$ and (d-f) $W=10$.
The initial state is a square distribution centered at $E_{\rm{mid}}=700$ and
of width $\Delta E=100$. Black lines are the moving time averages on the
logarithmic timescale with window $[\log_{10}t-\Delta t,\log_{10}t+\Delta t]$
and $\Delta t=0.02$. Green horizontal dashed line indicates the infinite time
average of the survival probability $S_{\infty}$. For $N=3$ and $4$, red lines
represent the analytic solution. For $N=4$, the grey points represent the
unaveraged value for comparison. The parameters used for each row correspond
to those that give the minimum kurtosis for $4$ particles in Fig.5, namely
$\gamma=30$ and $\tau=15$ for $W=2$, and $\gamma=20$ and $\tau=7.5$ for
$W=10$.

We take an initial state that has a homogeneous probability distribution
centered at $E_{\rm{mid}}$ inside an energy window of width $\Delta E$

$\rho(E)=\left\\{\\!\\!\\!\\!\begin{array}[]{cl}\frac{1}{2\Delta
E}&\mbox{for}\ E\in[E_{\rm{mid}}-\Delta E,E_{\rm{mid}}+\Delta E]\\\
0&\mbox{otherwise.}\end{array}\right.$

In [106, 95, 67] a general analytic solution was derived for the survival
probability for chaotic systems. For the square distribution used here, this
solution is given by

$S_{P}(t)=\frac{1-S_{\infty}}{\eta-1}\left[\eta\frac{\sin^{2}(\Delta
E\,t)}{(\Delta E\,t)^{2}}-b_{2}\left(\frac{\Delta
E\,t}{\pi\eta}\right)\right]+S_{\infty}\;.$ (8)

Here, $\eta$ is the number of energy eigenvectors in the energy window $\Delta
E$. The initial decay of the survival probability is captured by the first
term in the brackets in Eq. (8), after which the dynamics is described by the
two-level form factor

$\displaystyle b_{2}(\bar{t})$ $\displaystyle=$
$\displaystyle\left[1-2\bar{t}+\bar{t}\ln(2\bar{t}+1)\right]\Theta(1-\bar{t})$
$\displaystyle{}+\left[\bar{t}\ln\left(\frac{2\bar{t}+1}{2\bar{t}-1}\right)-1\right]\Theta(\bar{t}-1),$

with $\Theta(\bar{t})$ the Heaviside step function.

In Fig. 7 the survival probability is shown for different number of particles
and number of wells (grey points). The moving time average (black lines)
smooths the data and allows us to identify the correlation hole with its ramp
towards $S_{\infty}$. For $W=2$ and $N=2,3$ [Figs. 7 (a,b)] there is no
obvious indication of a dip below $S_{\infty}$ that could be described by the
two-level form factor. For $N=4$ [Fig. 7 (c)], a noticeable dip manifests
below the saturation point and follows closely the analytic solution in Eq.
(8) (red lines). This is indicative of chaotic behavior. In the $W=10$ well
system the correlation hole for $N=4$ is even more pronounced [Fig. 7 (f)]
with the initial decay and subsequent ramp of $S_{P}(t)$ matching precisely
the analytics. For $N=3$ and $W=10$ [Fig. 7 (e)], partial revivals still
obscure the minimum of the correlation hole, but the ramp towards saturation
can be seen to follow the analytic results.

###  4.4 Dependence on the number of wells

For a more systematic analysis of the onset of chaos as the number of wells in
the system is increased, we show in Fig. 8 (a) the minimum of the kurtosis
$\mathcal{K}_{\rm{min}}=\text{min}[\mathcal{K}_{T}]$ in the range
$\gamma,\tau\in[0,100]$. Here we fix the size of the box $L$ and focus on
$N=3$ and $N=4$ particles. For both cases the number of wells dictates the
emergence of chaos, but to different degrees. For the trivial case of $W=1$
(no barriers) the kurtosis is large and both $N=3$ and $N=4$ are integrable.
For $W>1$ the kurtosis of $N=4$ takes low values, $\mathcal{K}_{\rm{min}}\sim
3$, indicating chaos, essentially irrespective of the number of wells when
$W\lesssim 13$. However, for $N=3$ the inclusion of more barriers causes a
more subtle change to the kurtosis, which decreases slowly as more barriers
are introduced, attaining a minimum of $\mathcal{K}_{\rm{min}}\sim 4$ in the
region of $6\lesssim W\lesssim 13$. It is in this region that the $N=3$ system
displays chaotic signatures as discussed in the previous sections.
Interestingly, increasing the number of wells further ($W>13$) results in an
increase of the kurtosis and indications of chaos are diminished. A similar
tendency is seen for $N=4$, albeit in a less drastic manner, as the kurtosis
increases at a lower rate. For both $N=3$ and $N=4$ we expect that as the
limit of $N/W\rightarrow 0$ is approached the contest between the ordering of
the particles in the wells and their interactions is reduced and that the
system slowly returns to what we would see in the continuum: particles in a
box [90].

In Fig. 8 (b) and Fig. 8 (c) we show the optimal interactions and barrier
heights for achieving the minimum kurtosis shown in Fig. 8 (a). Figure 8 (c)
shows that the optimal barrier strengths $\tau$ for $N=3$ and $N=4$ have close
agreement for all $W$, and that these values decrease with increasing number
of wells (for $W\leq 20$). This reduction in $\tau$ is necessary to preserve
the competition between barriers and interactions, as when the number of
barriers is increased the impact of $\tau$ is magnified. This effect can be
seen in the shift of the chaotic region in Fig. 5 (c) and (f). Fig. 8 (b)
shows that the optimal interaction strengths $\gamma$ found for both $N=3$ and
$N=4$ converge to a similar value in the region $7\lesssim W\lesssim 20$,
which encompasses the areas of low kurtosis in Fig. 8 (a). This suggests that
both $N=3$ and $N=4$ particles are chaotic in the same regions of the
parameter space $[\tau,\gamma]$. Increasing the number of wells beyond $W=20$
breaks this trend and the parameters for $N=3$ and $N=4$ diverge as the
indications of chaos are lost.

Figure 8: (a) Minimum of the kurtosis as a function of the number of wells for
$N=3$ (red triangles) and $N=4$ (yellow circles). Kurtosis is calculated in an
energy window of width $\Delta E=100$ around $E_{\rm{mid}}=700$. (b)
Interaction and (c) barrier height corresponding to the points of minimum
kurtosis in (a).

##  5 Discussion and Conclusion

To summarize, between $N=2$ and $N=4$ our model makes a clear transition. For
$N=4$ particles and just $W=2$ wells, there are clear signatures of chaos when
both the barrier and the interactions have finite strength. The density of
states grows rapidly with energy, and numerical analysis gives evidence of the
highly-correlated spectrum typical of random matrices, of the validity of the
off-diagonal ETH, and of a clear ramp in the survival probability. As the
number of wells is increased, evidence for the onset of chaos becomes more
robust, even for weaker barriers and interactions.

In contrast, for $N=2$, the spectrum is numerically indistinguishable from an
integrable system throughout the parameter window for $W=2$, and deviates only
slightly from this as the number of wells is increased. Several possible
origins for the non-chaotic behavior of $N=2$ can be hypothesized, including
some undiscovered integral of motion, some sort of partial integrability of a
subset of eigenstates of the Hamiltonian, a proliferation of accidental
degeneracies from several sources that create the appearance of integrability,
or a combination of all these. These possibilities will be considered in more
detail in a subsequent work. For now, we note with some irony that
integrability is much less generic and more complicated than chaos.

The intermediate case of $N=3$ stands at the ragged edge where symmetry and
integrability dissolve into chaos and randomness. By tuning the interaction
and barrier parameters and increasing the number of barriers, the full gamut
of possibilities can be realized on the same small atomic system. The $N=3$
case has some common aspects with the $N=2$ case and other features similar to
$N=4$. For example, the density of states grows with the energy like the $N=4$
and unlike the $N=2$ case, although it grows sublinearly while the $N=4$ case
grows linearly. On the other hand, the probability distributions for the off-
diagonal elements of the kinetic energy are closer to the $N=2$ than to the
$N=4$ case for only two wells, but when increasing the number of wells, it
gets closer to the Gaussian distribution characteristic of the $N=4$ case. The
variety of possible scenarios for $N=3$ atoms is most clear in Fig. 8 (a),
where the degree of chaoticity measured with the minimum kurtosis
$\mathcal{K_{\mathrm{min}}}$ makes sweeping changes as the number of wells is
increased.

As we discuss above, current experiments with subwavelength lattice potentials
and deterministically prepared few-body states provide the ideal platform to
probe this boundary between integrability and chaos. Furthermore, these small
atomic systems have been proposed as the working units of larger quantum
information processing devices and protocols. Therefore, understanding their
information, control, and entanglement properties in these different regimes
becomes important.

For example, this model lends itself naturally to “digitization”. The number
of particles in each well becomes a useful observable in the infinite barrier
limit; they are integrals of motion in fact. Similarly, coherent
superpositions of eigenstates of these or other integrals of motion in the
limiting edge models could be used for storing quantum information (c.f.
[83]). A quench from an integrable limit to the chaotic parameter regime would
break these integrals of motion and effectively scramble the information held
in the initial state.

An important extension of this work that is experimentally relevant is
determining how sensitive our results are to the idealizations of delta-
barriers, precisely symmetric positions and uniform barrier height. We expect
that small deviations of the periodic Kronig-Penney lattice, such as finite
width barriers, would not significantly alter the chaotic regions of our
system. The versatility of our model also allows us to explore the possibility
of emergent integrability [1] when more disorder is introduced into the system
via non-regularly spaced barriers or barriers of different heights. Both
aspects are inspiring and we leave them for future research. Another
interesting scenario which we have not explored is the case in which the
interactions are attractive, whereby the system is now furnished with a bound
state and whose limit at infinite interactions is the so-called super Tonks-
Girardeau gas [3, 5, 45].

The authors thank M. Olshanii, T. Busch, A. Fabra and M. Boubakour for
insights on integrability and conversations about chaos. TF acknowledges
support from JSPS KAKENHI-21K13856 and the Okinawa Institute of Science and
Technology Graduate University. We are grateful for the help and support
provided by the Scientific Computing and Data Analysis section of Research
Support Division at OIST. LFS was supported by the NSF grant No. DMR-1936006.
M.A.G.M. acknowledges funding from the Spanish Ministry of Education and
Vocational Training (MEFP) through the Beatriz Galindo program 2018
(BEAGAL18/00203) and Spanish Ministry MINECO (FIDEUA
PID2019-106901GBI00/10.13039/501100011033).

##  Appendix A Density of states derivation

To derive the density of states, we first calculate the total number of states
(with any symmetry or parity) and energy less than $E$:

$\displaystyle\mathcal{N}(E)$ $\displaystyle=$
$\displaystyle\frac{1}{2^{N}}\frac{\pi^{N/2}}{\Gamma(N/2+1)}E^{N/2}$
$\displaystyle+O\left[N,W,\tau,\gamma\right](E^{(N-1)/2}).$

The density of states in Eq. (2) is the derivative of this with respect to
energy.

To establish this result (A), first consider the simplest limiting case
$H(N,W,0,0)$. There is a solution of $H(N,W,0,0)$ for every set of non-
negative integers ${\bf n}=\\{n_{1},\ldots,n_{N}\\}$ with energy
$E_{\mathbf{n}}=\sum_{i=1}^{N}n_{i}^{2}$. The space of solutions therefore is
a (hyper)cubic lattice in the all-positive ‘quadrant’ (really $2^{N}$-rant) of
$\mathbb{R}^{N}$. Since each state takes up a unit volume in this quantum
number space, to find the number of states $\mathcal{N}_{N}(E)$ with energy
less than $E$, one takes the volume of an $N$-ball with radius $r=\sqrt{E}$
and divides by $2^{N}$ to account for the all-positive condition, giving the
leading term in Eq. (A). The first correction term comes from the $N$-sphere
boundary of the $N$-ball, which has one dimension lower.

The spectrum of $H(2,2,0,0)$ is depicted as model 1 in Fig. 9. In this simple
case, states with energy less than $E$ lie within the quarter circle with
radius $r=\sqrt{E}$. That quarter disk therefore has area $\pi E/4$, agreeing
with (A).

Figure 9: These four diagrams depict the spectrum for the four solvable cases
of $N=2$ particle and $W=2$ wells; c.f. Fig. 1 and Sec. 3.1. Model 1: no
interactions or barriers; model 2: infinite interactions no barriers; model 3
no interactions, infinite barriers; model 4 infinite interactions and
barriers. Each dot represents an energy eigenstate with energy
$E=n_{1}^{2}+n_{2}^{2}$ given by the sum of the squares of the integer
coordinates $(n_{1},n_{2})$ of the point. States with $n_{1}\geq n_{2}$
represent the symmetrized states with positive parity in black and with with
negative parity in cyan. States with $n_{1}<n_{2}$ represented by magenta and
yellow dots are antisymmetric and have negative and positive parity,
respectively. In models 1 and 2, all states are two-fold degenerate except in
model 1 (no barriers, no interactions) when $n_{1}=n_{2}$. Additionally,
models 3 and 4 have additional two-fold degeneracies (in model 4, states along
the diagonal represented by split dots with twice the area), four-fold (in
model 4, represented by quartered dots along the diagonal) and eight-fold (in
models 3 and 4, represented by pairs of quartered dots with exchanged integer
quantum numbers). The dashed quarter-circle is the boundary $E=150$ with
radius $\sqrt{150}$ and is included to aid visualization of the spectral flow.

Parallel arguments give the same leading terms for the other three solvable
corner models. For example, in the limit of no barriers and infinite
interactions $H(N,W,0,\infty)$, we follow the construction of Girardeau and
find that there are $N!$ states for every strictly increasing set of integers
$\\{{\bf n}|n_{1}<n_{2}<\cdots n_{N}\\}$; c.f. [47]. Therefore, infinite
interactions exclude cases when two or more quantum numbers are the same. This
is depicted for $H(2,2,0,\infty)$ in Fig. 9, where the spectrum is missing the
diagonal states with $n_{1}=n_{2}$. For $N=2$, the number of states missing
from the area estimation in Eq. (A) is therefore proportional to the length of
this boundary $r=\sqrt{E}$. More generally for $N$ particles, the geometrical
structures in the quantum number ‘quadrant’ corresponding to all numbers
different remains $N$-dimensional even in the presence of infinite
interactions. However, the structures with two quantum numbers equal
correspond to interior boundaries of the quadrant and have dimension $N-1$.
Therefore corrections to account for the ‘missing states’ appear at the
subleading order $r^{N-1}=E^{(N-1)/2}$. Further corrections to the number of
states appear at next-to-subleading order $E^{(N-2)/2}$ when either two pairs
of quantum numbers are the same or three quantum numbers are the same.

Note that in the example with $N=2$ in Fig. 9, the Tonks-Girardeau map shifts
the symmetric states $(n_{1},n_{2})$ in model 1 to $(n_{1}+1,n_{2})$ in model
2. Since an integrable model connects these two limiting cases, this
establishes a one-to-one adiabatic mapping between the two spectra. From this
mapping the number of level crossings that occur as $\gamma$ is tuned from $0$
to $\infty$ can be explicitly calculated without actually solving for the
spectrum on the integrable model that links these two cases.

Similarly, for the case $H(N,W,\infty,0)$ of infinite barriers and no
interactions, there are $W^{N}$ solutions for every set of non-negative
integers ${\bf n}$ where all integers $n_{i}$ are multiples of the number of
wells $W$. This increases the volume associated with each set of quantum
numbers from $1$ to $W^{N}$ and that factor of $1/W^{N}$ exactly cancels the
degeneracy factor $W^{N}$ giving the same leading term. This is depicted for
the simplest case of $N=2$ and $W=2$ in model 3 of Fig. 9. As before, we do
expect the coefficient on the subleading term to depend on the intricate
combinatorics of putting $N$ identical particles in $W$ wells. A parallel
argument holds for the fourth solvable model $H(N,W,\infty,\infty)$.

Note again that these four models with exact solutions are connected by
integral models and exact spectral maps can be constructed for all of these
cases. Because these four extreme cases all have the same leading term in the
density of states that depends only on $N$, we assume that all the models that
lie in this region of parameter space have the same property, and this is
further supported by the phase space argument presented in Sec. 4.1.

When we elect to consider only one symmetry sector, e.g. bosons with positive
parity, that reduces the number of states (A) by a factor of $1/2$ for parity
and $1/N!$ for symmetrization at leading order. At subleading order, the
correction coefficient depends on the combinatorics of $N$ particles in $W$
wells and on the parameters $(\tau,\gamma)$. For example, for the simplest
case of $N=2$ and $W=2$ depicted, we see the importance of states along the
diagonal $n_{1}=n_{2}$, which would result in corrections of order $\sqrt{E}$
for the length of that diagonal.

As a final comment, if the leading term of the density of states is
independent of the parameters $(\tau,\gamma)$, then the subleading term
contains information about the density of level crossings for the integrable
models which presumably becomes level repulsions in the non-integrable
parameter region. Future work will investigate this connection between
spectral flow, level crossings, and energy level statistics.

## References

  * Abanin et al. [2019] Dmitry A. Abanin, Ehud Altman, Immanuel Bloch, and Maksym Serbyn.  Colloquium: Many-body localization, thermalization, and entanglement.  _Rev. Mod. Phys._ , 91:021001, May 2019.  doi: 10.1103/RevModPhys.91.021001.  URL https://link.aps.org/doi/10.1103/RevModPhys.91.021001.
  * Alhassid and Levine [1992] Y. Alhassid and R. D. Levine.  Spectral autocorrelation function in the statistical theory of energy levels.  _Phys. Rev. A_ , 46:4650–4653, 1992.  doi: 10.1103/PhysRevA.46.4650.  URL https://link.aps.org/doi/10.1103/PhysRevA.46.4650.
  * Astrakharchik et al. [2005] G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini.  Beyond the Tonks-Girardeau gas: Strongly correlated regime in quasi-one-dimensional bose gases.  _Phys. Rev. Lett._ , 95:190407, Nov 2005.  doi: 10.1103/PhysRevLett.95.190407.  URL https://link.aps.org/doi/10.1103/PhysRevLett.95.190407.
  * Barredo et al. [2016] Daniel Barredo, Sylvain de Léséleuc, Vincent Lienhard, Thierry Lahaye, and Antoine Browaeys.  An atom-by-atom assembler of defect-free arbitrary two-dimensional atomic arrays.  _Science_ , 354(6315):1021–1023, 2016.  doi: 10.1126/science.aah3778.
  * Batchelor et al. [2005] M T Batchelor, M Bortz, X W Guan, and N Oelkers.  Evidence for the super Tonks-Girardeau gas.  _J. Stat. Mech._ , 2005(10):L10001–L10001, oct 2005.  doi: 10.1088/1742-5468/2005/10/l10001.  URL https://doi.org/10.1088/1742-5468/2005/10/l10001.
  * Becher et al. [2020] J. H. Becher, E. Sindici, R. Klemt, S. Jochim, A. J. Daley, and P. M. Preiss.  Measurement of identical particle entanglement and the influence of antisymmetrization.  _Phys. Rev. Lett._ , 125:180402, Oct 2020.  doi: 10.1103/PhysRevLett.125.180402.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.180402.
  * Bergschneider et al. [2019] Andrea Bergschneider, Vincent M Klinkhamer, Jan Hendrik Becher, Ralf Klemt, Lukas Palm, Gerhard Zürn, Selim Jochim, and Philipp M Preiss.  Experimental characterization of two-particle entanglement through position and momentum correlations.  _Nat. Phys._ , 15(7):640–644, 2019.  doi: 10.1038/s41567-019-0508-6.
  * Bernien et al. [2017] Hannes Bernien, Sylvain Schwartz, Alexander Keesling, Harry Levine, Ahmed Omran, Hannes Pichler, Soonwon Choi, Alexander S. Zibrov, Manuel Endres, Markus Greiner, Vladan Vuletić, and Mikhail D. Lukin.  Probing many-body dynamics on a 51-atom quantum simulator.  _Nature_ , 551(7682):579–584, 2017.  ISSN 1476-4687.  doi: 10.1038/nature24622.  URL https://doi.org/10.1038/nature24622.
  * Berry [1981] M. V. Berry.  Quantizing a classically ergodic system: Sinai’s billiard and the KKR method.  _Ann. Phys._ , 131(1):163–216, January 1981.  ISSN 0003-4916.  doi: 10.1016/0003-4916(81)90189-5.  URL https://www.sciencedirect.com/science/article/pii/0003491681901895.
  * Berry and Tabor [1977] M. V. Berry and M. Tabor.  Level clustering in the regular spectrum.  _Proc. R. Soc. Lond. A_ , 356:375 – 394, 1977.  doi: 10.1098/rspa.1977.0140.
  * Berry and Wilkinson [1984] M. V. Berry and M. Wilkinson.  Diabolical Points in the Spectra of Triangles.  _Proc. R. Soc. Lond. A_ , 392(1802):15–43, 1984\.  ISSN 0080-4630.  doi: 10.1098/rspa.1984.0022.  URL https://www.jstor.org/stable/2397740.
  * Bertini et al. [2021] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič.  Finite-temperature transport in one-dimensional quantum lattice models.  _Rev. Mod. Phys._ , 93:025003, May 2021.  doi: 10.1103/RevModPhys.93.025003.  URL https://link.aps.org/doi/10.1103/RevModPhys.93.025003.
  * Beugeling et al. [2015] Wouter Beugeling, Roderich Moessner, and Masudul Haque.  Off-diagonal matrix elements of local operators in many-body quantum systems.  _Phys. Rev. E_ , 91:012144, Jan 2015.  doi: 10.1103/PhysRevE.91.012144.  URL https://link.aps.org/doi/10.1103/PhysRevE.91.012144.
  * Bloch et al. [2008] Immanuel Bloch, Jean Dalibard, and Wilhelm Zwerger.  Many-body physics with ultracold gases.  _Rev. Mod. Phys._ , 80:885–964, Jul 2008.  doi: 10.1103/RevModPhys.80.885.  URL https://link.aps.org/doi/10.1103/RevModPhys.80.885.
  * Blume [2012] D Blume.  Few-body physics with ultracold atomic and molecular systems in traps.  _Rep. Prog. Phys._ , 75(4):046401, mar 2012.  doi: 10.1088/0034-4885/75/4/046401.  URL https://doi.org/10.1088/0034-4885/75/4/046401.
  * Blume [2010] Doerte Blume.  Jumping from two and three particles to infinitely many.  _Physics_ , 3:74, 2010.  doi: 10.1103/Physics.3.74.
  * Borgonovi et al. [2016] F. Borgonovi, F. M. Izrailev, L. F. Santos, and V. G. Zelevinsky.  Quantum chaos and thermalization in isolated systems of interacting particles.  _Phys. Rep._ , 626:1, 2016.  doi: 10.1016/j.physrep.2016.02.005.
  * Bourgain et al. [2013] R. Bourgain, J. Pellegrino, A. Fuhrmanek, Y. R. P. Sortais, and A. Browaeys.  Evaporative cooling of a small number of atoms in a single-beam microscopic dipole trap.  _Phys. Rev. A_ , 88:023428, Aug 2013.  doi: 10.1103/PhysRevA.88.023428.  URL https://link.aps.org/doi/10.1103/PhysRevA.88.023428.
  * Brenes et al. [2020a] Marlon Brenes, John Goold, and Marcos Rigol.  Low-frequency behavior of off-diagonal matrix elements in the integrable XXZ chain and in a locally perturbed quantum-chaotic XXZ chain.  _Phys. Rev. B_ , 102:075127, Aug 2020a.  doi: 10.1103/PhysRevB.102.075127.  URL https://link.aps.org/doi/10.1103/PhysRevB.102.075127.
  * Brenes et al. [2020b] Marlon Brenes, Tyler LeBlond, John Goold, and Marcos Rigol.  Eigenstate thermalization in a locally perturbed integrable system.  _Phys. Rev. Lett._ , 125:070605, Aug 2020b.  doi: https://doi.org/10.1103/PhysRevLett.125.070605.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.070605.
  * Brenes et al. [2021] Marlon Brenes, Silvia Pappalardi, Mark T. Mitchison, John Goold, and Alessandro Silva.  Out-of-time-order correlations and the fine structure of eigenstate thermalisation.  _arXiv:2103.01161_ , 2021.  URL https://arxiv.org/abs/2103.01161.
  * Brody et al. [1981] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong.  Random-matrix physics: spectrum and strength fluctuations.  _Rev. Mod. Phys._ , 53:385, 1981.  doi: 10.1103/RevModPhys.53.385.  URL https://link.aps.org/doi/10.1103/RevModPhys.53.385.
  * Campbell et al. [2014] Steve Campbell, Miguel Ángel García-March, Thomás Fogarty, and Thomas Busch.  Quenching small quantum gases: Genesis of the orthogonality catastrophe.  _Phys. Rev. A_ , 90:013617, Jul 2014.  doi: 10.1103/PhysRevA.90.013617.  URL https://link.aps.org/doi/10.1103/PhysRevA.90.013617.
  * Cetina et al. [2016] Marko Cetina, Michael Jag, Rianne S. Lous, Isabella Fritsche, Jook T. M. Walraven, Rudolf Grimm, Jesper Levinsen, Meera M. Parish, Richard Schmidt, Michael Knap, and Eugene Demler.  Ultrafast many-body interferometry of impurities coupled to a Fermi sea.  _Science_ , 354(6308):96–99, 2016.  ISSN 0036-8075.  doi: 10.1126/science.aaf5134.
  * Cheinet et al. [2008] P. Cheinet, S. Trotzky, M. Feld, U. Schnorrberger, M. Moreno-Cardoner, S. Fölling, and I. Bloch.  Counting atoms using interaction blockade in an optical superlattice.  _Phys. Rev. Lett._ , 101:090404, Aug 2008.  doi: 10.1103/PhysRevLett.101.090404.  URL https://link.aps.org/doi/10.1103/PhysRevLett.101.090404.
  * Chenu et al. [2019] Aurélia Chenu, Javier Molina-Vilaplana, and Adolfo del Campo.  Work statistics, Loschmidt echo and information scrambling in chaotic quantum systems.  _Quantum_ , 3:127, March 2019.  ISSN 2521-327X.  doi: 10.22331/q-2019-03-04-127.  URL https://doi.org/10.22331/q-2019-03-04-127.
  * Cotler et al. [2017] J. Cotler, N. Hunter-Jones, J. Liu, and B. Yoshida.  Chaos, complexity, and random matrices.  _J. High Energy Phys._ , 2017(11):48, Nov 2017\.  doi: 10.1007/JHEP11(2017)048.  URL https://doi.org/10.1007/JHEP11(2017)048.
  * D’Alessio et al. [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol.  From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics.  _Adv. Phys._ , 65(3):239–362, 2016.  doi: 10.1080/00018732.2016.1198134.  URL https://doi.org/10.1080/00018732.2016.1198134.
  * de la Cruz et al. [2020] Javier de la Cruz, Sergio Lerma-Hernández, and Jorge G. Hirsch.  Quantum chaos in a system with high degree of symmetries.  _Phys. Rev. E_ , 102:032208, Sep 2020.  doi: 10.1103/PhysRevE.102.032208.  URL https://link.aps.org/doi/10.1103/PhysRevE.102.032208.
  * Endres et al. [2016] Manuel Endres, Hannes Bernien, Alexander Keesling, Harry Levine, Eric R. Anschuetz, Alexandre Krajenbrink, Crystal Senko, Vladan Vuletic, Markus Greiner, and Mikhail D. Lukin.  Atom-by-atom assembly of defect-free one-dimensional cold atom arrays.  _Science_ , 354(6315):1024–1027, 2016.  doi: 10.1126/science.aah3752.
  * Feldmeier and Schnack [2000] Hans Feldmeier and Jürgen Schnack.  Molecular dynamics for fermions.  _Rev. Mod. Phys._ , 72:655–688, Jul 2000.  doi: 10.1103/RevModPhys.72.655.  URL https://link.aps.org/doi/10.1103/RevModPhys.72.655.
  * Flambaum and Izrailev [1997] V. V. Flambaum and F. M. Izrailev.  Statistical theory of finite Fermi systems based on the structure of chaotic eigenstates.  _Phys. Rev. E_ , 56:5144, 1997.  doi: 10.1103/PhysRevE.56.5144.
  * Flambaum et al. [1994] V. V. Flambaum, A. A. Gribakina, G. F. Gribakin, and M. G. Kozlov.  Structure of compound states in the chaotic spectrum of the Ce atom: Localization properties, matrix elements, and enhancement of weak perturbations.  _Phys. Rev. A_ , 50:267–296, 1994.  doi: 10.1103/PhysRevA.50.267.
  * Flambaum et al. [1996] V. V. Flambaum, F. M. Izrailev, and G. Casati.  Towards a statistical theory of finite Fermi systems and compound states: Random two-body interaction approach.  _Phys. Rev. E_ , 54:2136–2139, Aug 1996.  doi: 10.1103/PhysRevE.54.2136.  URL https://link.aps.org/doi/10.1103/PhysRevE.54.2136.
  * Fogarty et al. [2020] Thomás Fogarty, Sebastian Deffner, Thomas Busch, and Steve Campbell.  Orthogonality catastrophe as a consequence of the quantum speed limit.  _Phys. Rev. Lett._ , 124:110601, Mar 2020.  doi: 10.1103/PhysRevLett.124.110601.  URL https://link.aps.org/doi/10.1103/PhysRevLett.124.110601.
  * García-March et al. [2014] M. A. García-March, B. Juliá-Díaz, G. E. Astrakharchik, J. Boronat, and A. Polls.  Distinguishability, degeneracy, and correlations in three harmonically trapped bosons in one dimension.  _Phys. Rev. A_ , 90:063605, Dec 2014.  doi: 10.1103/PhysRevA.90.063605.  URL https://link.aps.org/doi/10.1103/PhysRevA.90.063605.
  * García-March et al. [2015] M. A. García-March, A. Yuste, B. Juliá-Díaz, and A. Polls.  Mesoscopic superpositions of Tonks-Girardeau states and the Bose-Fermi mapping.  _Phys. Rev. A_ , 92:033621, Sep 2015.  doi: 10.1103/PhysRevA.92.033621.  URL https://link.aps.org/doi/10.1103/PhysRevA.92.033621.
  * Garcia-March et al. [2018] M A Garcia-March, S van Frank, M Bonneau, J Schmiedmayer, M Lewenstein, and Lea F Santos.  Relaxation, chaos, and thermalization in a three-mode model of a Bose–Einstein condensate.  _New J. Phys._ , 20(11):113039, nov 2018.  doi: 10.1088/1367-2630/aaed68.  URL https://doi.org/10.1088/1367-2630/aaed68.
  * Gaudin [1971] M. Gaudin.  Boundary energy of a Bose gas in one dimension.  _Physical Review A_ , 4(1):386–394, July 1971\.  doi: 10.1103/PhysRevA.4.386.  URL http://link.aps.org/doi/10.1103/PhysRevA.4.386.
  * Gaudin [2014] Michel Gaudin.  _The Bethe Wavefunction_.  Cambridge University Press, Trans. Caux, Jean-Sébastien edition, April 2014.  ISBN 978-1-107-04585-9.  doi: 10.1017/CBO9781107053885.
  * Girardeau [1960] Marvin Girardeau.  Relationship between Systems of Impenetrable Bosons and Fermions in One Dimension.  _Journal of Mathematical Physics_ , 1(6):516–523, November 1960.  ISSN 00222488.  doi: 10.1063/1.1703687.  URL http://jmp.aip.org/resource/1/jmapaq/v1/i6/p516_s1?isAuthorized=no.
  * Goold et al. [2011] J. Goold, T. Fogarty, N. Lo Gullo, M. Paternostro, and Th. Busch.  Orthogonality catastrophe as a consequence of qubit embedding in an ultracold Fermi gas.  _Phys. Rev. A_ , 84:063632, Dec 2011.  doi: 10.1103/PhysRevA.84.063632.  URL https://link.aps.org/doi/10.1103/PhysRevA.84.063632.
  * Guhr and Weidenmüller [1990] T. Guhr and H.A. Weidenmüller.  Correlations in anticrossing spectra and scattering theory. Analytical aspects.  _Chem. Phys._ , 146:21 – 38, 1990.  doi: 10.1016/0301-0104(90)90003-R.  URL http://www.sciencedirect.com/science/article/pii/030101049090003R.
  * Guhr et al. [1998] T. Guhr, A. Mueller-Gröeling, and H. A. Weidenmüller.  Random matrix theories in quantum physics: Common concepts.  _Phys. Rep._ , 299:189, 1998.  doi: 10.1016/S0370-1573(97)00088-4.
  * Haller et al. [2009] Elmar Haller, Mattias Gustavsson, Manfred J. Mark, Johann G. Danzl, Russell Hart, Guido Pupillo, and Hanns-Christoph Nägerl.  Realization of an excited, strongly correlated quantum gas phase.  _Science_ , 325(5945):1224–1227, 2009.  doi: 10.1126/science.1175850.
  * Harshman [2012] N. L. Harshman.  Symmetries of three harmonically trapped particles in one dimension.  _Phys. Rev. A_ , 86:052122, Nov 2012.  doi: 10.1103/PhysRevA.86.052122.  URL https://link.aps.org/doi/10.1103/PhysRevA.86.052122.
  * Harshman [2016a] N. L. Harshman.  One-dimensional traps, two-body interactions, few-body symmetries: II. $N$ particles.  _Few-Body Systems_ , 57(1):45–69, 2016a.  ISSN 0177-7963, 1432-5411.  doi: 10.1007/s00601-015-1025-5.  URL http://link.springer.com.proxyau.wrlc.org/article/10.1007/s00601-015-1025-5.
  * Harshman [2017] N. L. Harshman.  Infinite barriers and symmetries for a few trapped particles in one dimension.  _Phys. Rev. A_ , 95(5):053616, May 2017.  doi: 10.1103/PhysRevA.95.053616.  URL https://link.aps.org/doi/10.1103/PhysRevA.95.053616.
  * Harshman et al. [2017] N. L. Harshman, Maxim Olshanii, A. S. Dehkharghani, A. G. Volosniev, Steven Glenn Jackson, and N. T. Zinner.  Integrable families of hard-core particles with unequal masses in a one-dimensional harmonic trap.  _Phys. Rev. X_ , 7:041001, Oct 2017.  doi: 10.1103/PhysRevX.7.041001.  URL https://link.aps.org/doi/10.1103/PhysRevX.7.041001.
  * Harshman [2016b] N.L. Harshman.  One-dimensional traps, two-body interactions, few-body symmetries: I. One, two, and three particles.  _Few-Body Systems_ , 57(1):11–43, 2016b.  doi: 10.1007/s00601-015-1024-6.
  * Hartmann et al. [1991] U. Hartmann, H.A. Weidenmüller, and T. Guhr.  Correlations in anticrossing spectra and scattering theory: Numerical simulations.  _Chem. Phys._ , 150(3):311 – 320, 1991.  ISSN 0301-0104.  doi: http://dx.doi.org/10.1016/0301-0104(91)87105-5.  URL http://www.sciencedirect.com/science/article/pii/0301010491871055.
  * He et al. [2010] Xiaodong He, Peng Xu, Jin Wang, and Mingsheng Zhan.  High efficient loading of two atoms into a microscopic optical trap by dynamically reshaping the trap with a spatial light modulator.  _Opt. Express_ , 18(13):13586–13592, 2010.  doi: 10.1364/OE.18.013586.
  * Heilmann and Lieb [1971] Ole J. Heilmann and Elliott H. Lieb.  Violation of the noncrossing rule: The Hubbard Hamiltonian for benzene.  _Ann. N.Y. Acad. Sci._ , 172(15):584–617, 1971\.  ISSN 1749-6632.  doi: 10.1111/j.1749-6632.1971.tb34956.x.  URL https://nyaspubs.onlinelibrary.wiley.com/doi/abs/10.1111/j.1749-6632.1971.tb34956.x.
  * Huber et al. [2021] David Huber, Oleksandr V. Marchukov, Hans-Werner Hammer, and Artem G. Volosniev.  Morphology of three-body quantum states from machine learning.  _New Journal of Physics_ , 2021.  doi: 10.1088/1367-2630/ac0576.  URL http://iopscience.iop.org/article/10.1088/1367-2630/ac0576.
  * Žnidarič [2020] Marko Žnidarič.  Weak integrability breaking: Chaos with integrability signature in coherent diffusion.  _Phys. Rev. Lett._ , 125:180605, Oct 2020.  doi: 10.1103/PhysRevLett.125.180605.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.180605.
  * Ivrii [2016] Victor Ivrii.  100 years of Weyl’s law.  _Bull. Math. Sci._ , 6(3):379–452, October 2016\.  ISSN 1664-3607, 1664-3615.  doi: 10.1007/s13373-016-0089-y.  URL https://link.springer.com/article/10.1007/s13373-016-0089-y.
  * Izrailev [2001] F. M. Izrailev.  Quantum-classical correspondence for isolated systems of interacting particles: Localization and ergodicity in energy space.  _Phys. Scr._ , T90:95–104, 2001.  doi: 10.1238/Physica.Topical.090a00095.
  * Izrailev [1990] Felix M. Izrailev.  Simple models of quantum chaos: Spectrum and eigenfunctions.  _Phys. Rep._ , 196(5):299 – 392, 1990.  ISSN 0370-1573.  doi: https://doi.org/10.1016/0370-1573(90)90067-C.  URL http://www.sciencedirect.com/science/article/pii/037015739090067C.
  * Jain [2008] Sudhir R. Jain.  Exact solution of the Schrödinger equation for a particle in a regular N-simplex.  _Phys. Lett. A_ , 372(12):1978–1981, March 2008\.  ISSN 0375-9601.  doi: 10.1016/j.physleta.2007.11.016.  URL https://www.sciencedirect.com/science/article/pii/S0375960107016027.
  * Kaufman et al. [2012] A. M. Kaufman, B. J. Lester, and C. A. Regal.  Cooling a single atom in an optical tweezer to its quantum ground state.  _Phys. Rev. X_ , 2:041014, Nov 2012.  doi: 10.1103/PhysRevX.2.041014.
  * Kaufman et al. [2014] A. M. Kaufman, B. J. Lester, C. M. Reynolds, M. L. Wall, M. Foss-Feig, K. R. A. Hazzard, A. M. Rey, and C. A. Regal.  Two-particle quantum interference in tunnel-coupled optical tweezers.  _Science_ , 345(6194):306–309, 2014.  doi: 10.1126/science.1250057.
  * Keller and Fogarty [2016] Tim Keller and Thomás Fogarty.  Probing the out-of-equilibrium dynamics of two interacting atoms.  _Phys. Rev. A_ , 94:063620, Dec 2016.  doi: 10.1103/PhysRevA.94.063620.  URL https://link.aps.org/doi/10.1103/PhysRevA.94.063620.
  * Knap et al. [2012] Michael Knap, Aditya Shashi, Yusuke Nishida, Adilet Imambekov, Dmitry A. Abanin, and Eugene Demler.  Time-dependent impurity in ultracold fermions: Orthogonality catastrophe and beyond.  _Phys. Rev. X_ , 2:041020, Dec 2012.  doi: 10.1103/PhysRevX.2.041020.  URL https://link.aps.org/doi/10.1103/PhysRevX.2.041020.
  * Krishnamurthy et al. [1982] H. R. Krishnamurthy, H. S. Mani, and H. C. Verma.  Exact solution of the schrodinger equation for a particle in a tetrahedral box.  _J. Phys. A_ , 15(7):2131–2137, July 1982.  ISSN 0305-4470.  doi: 10.1088/0305-4470/15/7/024.  URL https://doi.org/10.1088/0305-4470/15/7/024.
  * Łacki et al. [2016] M. Łacki, M. A. Baranov, H. Pichler, and P. Zoller.  Nanoscale “dark state” optical potentials for cold atoms.  _Phys. Rev. Lett._ , 117:233001, Nov 2016.  doi: 10.1103/PhysRevLett.117.233001.  URL https://link.aps.org/doi/10.1103/PhysRevLett.117.233001.
  * LeBlond et al. [2019] Tyler LeBlond, Krishnanand Mallayya, Lev Vidmar, and Marcos Rigol.  Entanglement and matrix elements of observables in interacting integrable systems.  _Phys. Rev. E_ , 100:062134, Dec 2019.  doi: 10.1103/PhysRevE.100.062134.  URL https://link.aps.org/doi/10.1103/PhysRevE.100.062134.
  * Lerma-Hernández et al. [2019] S. Lerma-Hernández, D. Villaseñor, M. A. Bastarrachea-Magnani, E. J. Torres-Herrera, L. F. Santos, and J. G. Hirsch.  Dynamical signatures of quantum chaos and relaxation time scales in a spin-boson system.  _Phys. Rev. E_ , 100:012218, Jul 2019.  doi: 10.1103/PhysRevE.100.012218.  URL https://link.aps.org/doi/10.1103/PhysRevE.100.012218.
  * Lester et al. [2015] Brian J. Lester, Niclas Luick, Adam M. Kaufman, Collin M. Reynolds, and Cindy A. Regal.  Rapid production of uniformly filled arrays of neutral atoms.  _Phys. Rev. Lett._ , 115:073003, Aug 2015.  doi: 10.1103/PhysRevLett.115.073003.  URL https://link.aps.org/doi/10.1103/PhysRevLett.115.073003.
  * Leviandier et al. [1986] Luc Leviandier, Maurice Lombardi, Rémi Jost, and Jean Paul Pique.  Fourier transform: A tool to measure statistical level properties in very complex spectra.  _Phys. Rev. Lett._ , 56:2449–2452, Jun 1986.  doi: 10.1103/PhysRevLett.56.2449.
  * Lewenstein et al. [2012] Maciej Lewenstein, Anna Sanpera, and Veronica Ahufinger.  _Ultracold Atoms in Optical Lattices: Simulating quantum many-body systems_.  Oxford University Press, 2012.  doi: 10.1093/acprof:oso/9780199573127.001.0001.
  * Leyvraz et al. [2013] F Leyvraz, A García, H Kohler, and T H Seligman.  Fidelity under isospectral perturbations: a random matrix study.  _J. Phys. A_ , 46(27):275303, 2013.  doi: 10.1088/1751-8113/46/27/275303.  URL http://stacks.iop.org/1751-8121/46/i=27/a=275303.
  * Lieb and Liniger [1963] Elliott H. Lieb and Werner Liniger.  Exact analysis of an interacting Bose gas. I. The general solution and the ground state.  _Phys. Rev._ , 130:1605–1616, May 1963.  doi: 10.1103/PhysRev.130.1605.  URL https://link.aps.org/doi/10.1103/PhysRev.130.1605.
  * Lombardi and Seligman [1993] M. Lombardi and T. H. Seligman.  Universal and nonuniversal statistical properties of levels and intensities for chaotic Rydberg molecules.  _Phys. Rev. A_ , 47:3571–3586, May 1993.  doi: 10.1103/PhysRevA.47.3571.  URL http://link.aps.org/doi/10.1103/PhysRevA.47.3571.
  * Luitz and Lev [2017] D. Luitz and Y. Bar Lev.  The ergodic side of the many-body localization transition.  _Ann. Phys.(Berlin)_ , 529(7):1600350, 2017.  doi: 10.1002/andp.201600350.
  * Maldacena et al. [2016] Juan Maldacena, Stephen H. Shenker, and Douglas Stanford.  A bound on chaos.  _J. High Energy Phys._ , 2016(8):106, 2016.  doi: 10.1007/JHEP08(2016)106.  URL http://dx.doi.org/10.1007/JHEP08(2016)106.
  * Mehta [1991] M. L. Mehta.  _Random Matrices_.  Academic Press, Boston, USA, 1991.
  * Mirkin and Wisniacki [2021] Nicolás Mirkin and Diego Wisniacki.  Quantum chaos, equilibration, and control in extremely short spin chains.  _Phys. Rev. E_ , 103:L020201, Feb 2021.  doi: 10.1103/PhysRevE.103.L020201.  URL https://link.aps.org/doi/10.1103/PhysRevE.103.L020201.
  * Murmann et al. [2015] Simon Murmann, Andrea Bergschneider, Vincent M Klinkhamer, Gerhard Zürn, Thomas Lompe, and Selim Jochim.  Two fermions in a double well: Exploring a fundamental building block of the Hubbard model.  _Phys. Rev. Lett._ , 114(8):080402, 2015.  doi: 10.1103/PhysRevLett.114.080402.
  * Nandkishore and Huse [2015] R. Nandkishore and D.A. Huse.  Many-body localization and thermalization in quantum statistical mechanics.  _Annu. Rev. Condens. Matter Phys._ , 6:15, 2015.  doi: 10.1146/annurev-conmatphys-031214-014726.
  * Oelkers et al. [2006] N. Oelkers, M. T. Batchelor, M. Bortz, and X.-W. Guan.  Bethe ansatz study of one-dimensional Bose and Fermi gases with periodic and hard wall boundary conditions.  _J. Phys. A_ , 39(5):1073, 2006.  ISSN 0305-4470.  doi: 10.1088/0305-4470/39/5/005.  URL http://stacks.iop.org/0305-4470/39/i=5/a=005.
  * Olshanii [1998] M. Olshanii.  Atomic scattering in the presence of an external confinement and a gas of impenetrable bosons.  _Phys. Rev. Lett._ , 81:938–941, Aug 1998.  doi: 10.1103/PhysRevLett.81.938.  URL https://link.aps.org/doi/10.1103/PhysRevLett.81.938.
  * Olshanii and Jackson [2015] Maxim Olshanii and Steven G. Jackson.  An exactly solvable quantum four-body problem associated with the symmetries of an octacube.  _New J. Phys._ , 17(10):105005, 2015.  ISSN 1367-2630.  doi: 10.1088/1367-2630/17/10/105005.  URL http://stacks.iop.org/1367-2630/17/i=10/a=105005.
  * Olshanii et al. [2018] Maxim Olshanii, Thibault Scoquart, Dmitry Yampolsky, Vanja Dunjko, and Steven Glenn Jackson.  Creating entanglement using integrals of motion.  _Phys. Rev. A_ , 97(1):013630, January 2018.  doi: 10.1103/PhysRevA.97.013630.  URL https://link.aps.org/doi/10.1103/PhysRevA.97.013630.
  * Pagano et al. [2014] Guido Pagano, Marco Mancini, Giacomo Cappellini, Pietro Lombardi, Florian Schäfer, Hui Hu, Xia-Ji Liu, Jacopo Catani, Carlo Sias, Massimo Inguscio, and Leonardo Fallani.  A one-dimensional liquid of fermions with tunable spin.  _Nat. Phys._ , 10(3):198–201, 2014.  doi: 10.1038/nphys2878.  URL https://doi.org/10.1038/nphys2878.
  * Pandey and Ramaswamy [1991] Akhilesh Pandey and Ramakrishna Ramaswamy.  Level spacings for harmonic-oscillator systems.  _Phys. Rev. A_ , 43:4237–4243, Apr 1991.  doi: 10.1103/PhysRevA.43.4237.  URL http://link.aps.org/doi/10.1103/PhysRevA.43.4237.
  * Pique et al. [1987] J. P. Pique, Y. Chen, R. W. Field, and J. L. Kinsey.  Chaos and dynamics on 0.5–300 ps time scales in vibrationally excited acetylene: Fourier transform of stimulated-emission pumping spectrum.  _Phys. Rev. Lett._ , 58:475–478, Feb 1987.  doi: 10.1103/PhysRevLett.58.475.  URL http://link.aps.org/doi/10.1103/PhysRevLett.58.475.
  * Prange [1997] R. E. Prange.  The spectral form factor is not self-averaging.  _Phys. Rev. Lett._ , 78:2280, 1997.  doi: 10.1103/PhysRevLett.78.2280.  URL https://link.aps.org/doi/10.1103/PhysRevLett.78.2280.
  * Regal [2016] Cindy Regal.  Bringing order to neutral atom arrays.  _Science_ , 354(6315):972–973, 2016.  doi: 10.1126/science.aaj2145.
  * Reshodko et al. [2019] Irina Reshodko, Albert Benseny, Judit Romhányi, and Thomas Busch.  Topological states in the Kronig–Penney model with arbitrary scattering potentials.  _New J. Phys._ , 21(1):013010, jan 2019.  doi: 10.1088/1367-2630/aaf9bf.  URL https://doi.org/10.1088/1367-2630/aaf9bf.
  * Rigol and Muramatsu [2005] Marcos Rigol and Alejandro Muramatsu.  Ground-state properties of hard-core bosons confined on one-dimensional optical lattices.  _Phys. Rev. A_ , 72:013604, Jul 2005.  doi: 10.1103/PhysRevA.72.013604.  URL https://link.aps.org/doi/10.1103/PhysRevA.72.013604.
  * Roberts et al. [2015] Daniel A. Roberts, Douglas Stanford, and Leonard Susskind.  Localized shocks.  _JHEP_ , 2015(3):51, 2015.  ISSN 1029-8479.  doi: 10.1007/JHEP03(2015)051.  URL http://dx.doi.org/10.1007/JHEP03(2015)051.
  * Santos et al. [2020] Lea F. Santos, Francisco Pérez-Bernal, and E. Jonathan Torres-Herrera.  Speck of chaos.  _Phys. Rev. Research_ , 2:043034, Oct 2020.  doi: 10.1103/PhysRevResearch.2.043034.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.043034.
  * Saskin et al. [2019] S. Saskin, J. T. Wilson, B. Grinkemeyer, and J. D. Thompson.  Narrow-line cooling and imaging of ytterbium atoms in an optical tweezer array.  _Phys. Rev. Lett._ , 122:143002, Apr 2019.  doi: 10.1103/PhysRevLett.122.143002.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.143002.
  * Schiulaz et al. [2018] M. Schiulaz, M. Távora, and L. F. Santos.  From few- to many-body quantum systems.  _Quantum Sci. Technol._ , 3:044006, 2018.  doi: 10.1088/2058-9565/aad913.
  * Schiulaz et al. [2019] Mauro Schiulaz, E. Jonathan Torres-Herrera, and Lea F. Santos.  Thouless and relaxation time scales in many-body quantum systems.  _Phys. Rev. B_ , 99:174313, May 2019.  doi: 10.1103/PhysRevB.99.174313.  URL https://link.aps.org/doi/10.1103/PhysRevB.99.174313.
  * Schiulaz et al. [2020] Mauro Schiulaz, E. Jonathan Torres-Herrera, Francisco Pérez-Bernal, and Lea F. Santos.  Self-averaging in many-body quantum systems out of equilibrium: Chaotic systems.  _Phys. Rev. B_ , 101:174312, May 2020.  doi: 10.1103/PhysRevB.101.174312.  URL https://link.aps.org/doi/10.1103/PhysRevB.101.174312.
  * Schnack and Feldmeier [1996] J. Schnack and H. Feldmeier.  Statistical properties of fermionic molecular dynamics.  _Nuc. Phys. A_ , 601(2):181 – 194, 1996.  ISSN 0375-9474.  doi: 10.1016/0375-9474(95)00505-6.  URL http://www.sciencedirect.com/science/article/pii/0375947495005056.
  * Schreiber et al. [2015] M. Schreiber, S. S. Hodgman, Pr. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch.  Observation of many-body localization of interacting fermions in a quasirandom optical lattice.  _Science_ , 349(6250):842–845, 2015.  ISSN 0036-8075.  doi: 10.1126/science.aaa7432.  URL http://science.sciencemag.org/content/349/6250/842.
  * Schymik et al. [2020] Kai-Niklas Schymik, Vincent Lienhard, Daniel Barredo, Pascal Scholl, Hannah Williams, Antoine Browaeys, and Thierry Lahaye.  Enhanced atom-by-atom assembly of arbitrary tweezer arrays.  _Phys. Rev. A_ , 102:063107, Dec 2020.  doi: 10.1103/PhysRevA.102.063107.
  * Serwane et al. [2011] Friedhelm Serwane, Gerhard Zürn, Thomas Lompe, TB Ottenstein, AN Wenz, and S Jochim.  Deterministic preparation of a tunable few-fermion system.  _Science_ , 332(6027):336–338, 2011.  doi: 10.1126/science.1201351.
  * Shaw [1974] G. B. Shaw.  Degeneracy in the particle-in-a-box problem.  _J. Phys. A_ , 7(13):1537–1546, September 1974\.  ISSN 0301-0015.  doi: 10.1088/0305-4470/7/13/008.  URL https://doi.org/10.1088/0305-4470/7/13/008.
  * Sowiński and García-March [2019] Tomasz Sowiński and Miguel Ángel García-March.  One-dimensional mixtures of several ultracold atoms: a review.  _Rep. Progr. Phys._ , 82(10):104401, 2019.  doi: 10.1088/1361-6633/ab3a80.
  * Stöckmann [2006] H-J. Stöckmann.  _Quantum Chaos: an introduction_.  Cambridge University Press, Cambridge, UK, 2006.  doi: 10.1017/CBO9780511524622.
  * Torres-Herrera and Santos [2017a] E. J. Torres-Herrera and Lea F. Santos.  Extended nonergodic states in disordered many-body quantum systems.  _Ann. Phys. (Berlin)_ , 529:1600284, 2017a.  ISSN 1521-3889.  doi: 10.1002/andp.201600284.  URL http://dx.doi.org/10.1002/andp.201600284.
  * Torres-Herrera and Santos [2017b] E. J. Torres-Herrera and Lea F. Santos.  Dynamical manifestations of quantum chaos: correlation hole and bulge.  _Philos. Trans. Royal Soc. A_ , 375(2108):20160434, 2017b.  doi: 10.1098/rsta.2016.0434.
  * Torres-Herrera et al. [2018] E. J. Torres-Herrera, Antonio M. García-García, and Lea F. Santos.  Generic dynamical features of quenched interacting quantum systems: Survival probability, density imbalance, and out-of-time-ordered correlator.  _Phys. Rev. B_ , 97:060303, Feb 2018.  doi: 10.1103/PhysRevB.97.060303.  URL https://link.aps.org/doi/10.1103/PhysRevB.97.060303.
  * Turner et al. [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić.  Weak ergodicity breaking from quantum many-body scars.  _Nat. Phys._ , 14(7):745–749, 2018.  ISSN 1745-2481.  doi: 10.1038/s41567-018-0137-5.  URL https://doi.org/10.1038/s41567-018-0137-5.
  * Turner [1984] J. W. Turner.  On the quantum particle in a polyhedral box.  _J. Phys. A_ , 17(14):2791–2797, October 1984\.  ISSN 0305-4470.  doi: 10.1088/0305-4470/17/14/022.  URL https://doi.org/10.1088/0305-4470/17/14/022.
  * Wang et al. [2018] Y. Wang, S. Subhankar, P. Bienias, M. Łacki, T-C. Tsui, M. A. Baranov, A. V. Gorshkov, P. Zoller, J. V. Porto, and S. L. Rolston.  Dark state optical lattice with a subwavelength spatial structure.  _Phys. Rev. Lett._ , 120:083601, Feb 2018.  doi: 10.1103/PhysRevLett.120.083601.  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.083601.
  * Wenz et al. [2013] AN Wenz, G Zürn, Simon Murmann, I Brouzos, T Lompe, and S Jochim.  From few to many: Observing the formation of a Fermi sea one atom at a time.  _Science_ , 342(6157):457–460, 2013.  doi: 10.1126/science.1240516.
  * Wilkie and Brumer [1991] Joshua Wilkie and Paul Brumer.  Time-dependent manifestations of quantum chaos.  _Phys. Rev. Lett._ , 67:1185–1188, Sep 1991.  doi: 10.1103/PhysRevLett.67.1185.  URL https://link.aps.org/doi/10.1103/PhysRevLett.67.1185.
  * Will et al. [2010] Sebastian Will, Thorsten Best, Ulrich Schneider, Lucia Hackermüller, Dirk-Sören Lühmann, and Immanuel Bloch.  Time-resolved observation of coherent multi-body interactions in quantum phase revivals.  _Nature_ , 465(7295):197–201, May 2010.  ISSN 1476-4687.  doi: 10.1038/nature09036.  URL https://doi.org/10.1038/nature09036.
  * Zangara et al. [2013] Pablo R. Zangara, Axel D. Dente, E. J. Torres-Herrera, Horacio M. Pastawski, A. Iucci, and Lea F. Santos.  Time fluctuations in isolated quantum systems of interacting particles.  _Phys. Rev. E_ , 88:032913, 2013.  doi: 10.1103/PhysRevE.88.032913.
  * Zelevinsky et al. [1996] V. Zelevinsky, B. A. Brown, N. Frazier, and M. Horoi.  The nuclear shell model as a testing ground for many-body quantum chaos.  _Phys. Rep._ , 276:85–176, 1996.  doi: 10.1016/S0370-1573(96)00007-5.
  * Zisling et al. [2021] Guy Zisling, Lea F. Santos, and Yevgeny Bar Lev.  How many particles make up a chaotic many-body quantum system?  _SciPost Phys._ , 10:88, 2021.  doi: 10.21468/SciPostPhys.10.4.088.  URL https://scipost.org/10.21468/SciPostPhys.10.4.088.