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

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.


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) [17,28,114], hinders localization [74,79,98], facilitates the phenomenon of many-body quantum scarring [8,107], leads to diffusive transport [12,55] and causes the fast spread of quantum information [75,91].
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 meanfield 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 [31,32,34,57,97].
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.

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 manybody systems [15,16,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 struc-ture 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 [36,46,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 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 [4,30,99] 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 po-tential. In this experiment, a few 6 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 [7,78] and multi-well systems [6].

Model
In the simplest case of equally-spaced barriers, the Hamiltonian describing our system takes the form This model realizes a system with N identical interacting particles of mass m trapped in a onedimensional box of length L that is disrupted by W − 1 delta-barriers. In the equation above, x i ∈ [0, π] are the positions of the particles scaled by the length L/π, the energy scale is provided by 1 = 2 π 2 /(2mL 2 ) (henceforth set to unity), and τ ≥ 0 and γ ≥ 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, τ, γ) is to understand how the signatures of quantum chaos scale with the number of particles N , the number of wells W , the barrier strength τ and the interaction strength γ. The particularly simple form (1) of the Hamiltonian means that it is amenable to analytic and numeric calculations. 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. [59,64,82,108].

Solvable and integrable limiting cases
Beyond these four exactly solvable special cases, the Hamiltonian H(N, W, τ, γ) is also integrable along the four edges of the (τ, γ) parameter space window (cf. Fig. 1). Referring to the four 'corner' models denoted above, the integrable limits are:  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.

Symmetries and degeneracies
To analyze signatures of chaos, the Hilbert space 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 (τ, γ) parameter window. First, it is symmetric under particle permutations S N , i.e. any permutation p ∈ S N is represented as a linear transformation in configuration (1) invariant. The second symmetry is total parity inversion Π about the center of the well, implemented as the affine linear transformation These two symmetries allow eigenstates to be classified by irreducible representations of S N and of parity Π and reduce the total Hilbert space H into subspaces with a given symmetry (cf. [47]). In this work, we focus on the sector of Hilbert space H [N ]+ ⊂ 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 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 H [N ]+ upon which V N,W vanishes when N is even or both N and W are odd. Also, in the limits τ → ∞ and γ → ∞, 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 (τ, γ) are particle exchanges S N and parity Π, 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 τ or γ, 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 τ and γ 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 ≥ 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 [48,101], 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.

Density of states
The density of states of the Hamiltonian H(N, W, τ, γ) is used to interpret some of the numerical results below. Perhaps surprisingly, it is independent of τ and γ to leading order in the energy E. The density of states in the sector H [N ]+ is (2) for the respective particle numbers. Symbols denote different system parameters: γ = 10 and τ = 0 (circles), γ = τ = 10 and W = 2 (triangles) and γ = τ = 10 and W = 10 (stars).
where we have used the notation O [N, W, τ, γ] 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. 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 τ or interaction strength γ. 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 (τ, γ) 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.

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 H [N ] + [46]. In the following we discuss the different indicators of chaos that are used in this work along with the numerical results.

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 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,38,85] 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 [44,76], P WD (s) = (πs/2) exp −πs 2 /4 , 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 For weak interactions, γ = 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 [53,113]. 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 (γ ≥ 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 (τ, γ); 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]), P β (s) = (β + 1)bs β exp(−bs β+1 ), For the Wigner-Dyson distribution, β ∼ 1, while the Poissonian distribution leads to β ∼ 0.
Close to the integrability limits, vanishing interactions (γ ≈ 0) with finite barriers (τ = 0) or finite interactions (γ = 0) with vanishing barrier height (τ ≈ 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 β < 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 (γ 0) and the barrier heights (τ 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 β remaining below 0.35.
For more particles, the distributions are different, with indications of energy level repulsion emerging for N = 3 in the region 15 γ 45 and 5 τ 35 where β ∼ 0.5 [ Fig. 4 (b)]. For N = 4 the fitting parameter β in almost the entire parameter space is larger than 0.5, with the areas of integrability confined to the edges of the (τ, γ) parameter space [ Fig. 4 (c)]. In fact, we find the maximum to be around β ∼ 0.83 at γ = 20 and τ = 12.5, indicating that the energy spacing distribution gets close to Wigner-Dyson.
In Figs Comparing the case of W = 10 to W = 2 in Fig. 4, we see that larger values of β are indeed found for W = 10 and the region of the parameter space where β 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 β ≈ 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.

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 to reach thermal equilibrium. In the equation above, C m = m|Ψ(0) is the overlap between the eigenstate |m of the Hamiltonian that describes the system and the initial state |Ψ(0) , and O mn = m|Ô|n . 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 offdiagonal 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 infinitetime average 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, 19-21, 66, 92], 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, where . indicates the average over all pairs of eigenstates |m = |n and σ is the standard de-viation of the distribution. For Gaussian distributions the kurtosis is given by KÔ = 3.
In our calculations of KÔ, 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 ∈ [E mid − ∆E, E mid + ∆E], where the center of the energy window is high in the spectrum at E mid = 700 and its width is ∆E = 100. In Figs. 5 (a)-(f) we show the inverse of the kurtosis, 1/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 τ 30 and interactions 20 γ 60, which coincides with the region of β > 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 K T 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 K T ≈ 3.3 for N = 4 particles in W = 10 wells [ Fig. 5 (f)]. Indeed, for 1 τ < 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 τ , 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 To quantify how the off-diagonal elements T mn behave as a function of the energy difference ω = |E m − E n | we also show in Fig. 6 the ratio which is equal to π/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 Γ T for N = 4 sit very close to π/2 over a large range of ω, 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 π/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 Γ T ≈ π/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.

Survival Probability
Spectral correlations get manifested also in the evolution of the survival probability, in the form of what is known as the correlation hole [2,29,43,51,67,69,71,73,86,95,[104][105][106]111], recently referred to also as the "ramp" [27]. The correlation hole corresponds to a dip below S ∞ = n |C n | 4 , which is the infinite-time average (saturation value) of S P (t) . This dip emerges also in the spectral form factor n,m e −i(En−Em)t , but contrary to this one, the survival probability is a true dynamical quantity. In the equation above, . indicates averages. The survival probability is non-selfaveraging [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 [26,35,42,63], and can be experimentally measured using interferometric techniques [24].
We take an initial state that has a homogeneous probability distribution centered at E mid inside an energy window of width ∆E In [67,95,106] a general analytic solution was derived for the survival probability for chaotic systems. For the square distribution used here, this solution is given by Here, η is the number of energy eigenvectors in the energy window ∆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 b 2 (t) = 1 − 2t +t ln(2t + 1) Θ(1 −t) with Θ(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 ∞ . For W = 2 and N = 2, 3 [ Figs. 7 (a,b)] there is no obvious indication of a dip below S ∞ 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 de-cay 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.

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 K min = min[K T ] in the range γ, τ ∈ [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, K min ∼ 3, indicating chaos, essentially irrespective of the number of wells when W 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 intro-duced, attaining a minimum of K min ∼ 4 in the region of 6 W 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 → 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 τ for N = 3 and N = 4 have close agreement for all W , and that these values decrease with increasing number of wells (for W ≤ 20). This reduction in τ is necessary to preserve the competition between barriers and interactions, as when the number of barriers is increased the impact of τ 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 γ found for both N = 3 and N = 4 converge to a similar value in the region 7 W 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 [τ, γ]. 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.

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 K 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 deltabarriers, 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]

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: The density of states in Eq. (2) is the derivative of this with respect to energy. To establish this result (9), 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 nonnegative integers n = {n 1 , . . . , n N } with energy E n = N i=1 n 2 i . The space of solutions therefore is a (hyper)cubic lattice in the all-positive 'quadrant' (really 2 N -rant) of R N . Since each state takes up a unit volume in this quantum number space, to find the number of states N N (E) with energy less than E, one takes the volume of an N -ball with radius r = √ E and divides by 2 N to account for the all-positive condition, giving the leading term in Eq. (9). 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 = √ E. That quarter disk therefore has area πE/4, agreeing with (9). Each dot represents an energy eigenstate with energy E = n 2 1 + 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 ≥ 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 √ 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, ∞), we follow the construction of Girardeau and find that there are N ! states for every strictly increasing set of integers {n|n 1 < n 2 < · · · 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, ∞) 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. (9) is therefore proportional to the length of this boundary r = √ 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 γ is tuned from 0 to ∞ 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, ∞, 0) of infinite barriers and no interactions, there are W N solutions for every set of non-negative integers 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 combina-torics of putting N identical particles in W wells. A parallel argument holds for the fourth solvable model H (N, W, ∞, ∞).
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 (9) 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 (τ, γ). 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 √ 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 (τ, γ), 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.