# From integrability to chaos: the quantum-classical correspondence in a
triple well bosonic model

Erick R. Castro<EMAIL_ADDRESS>Centro Brasileiro de Pesquisas Físicas/MCTI,
22290-180, Rio de Janeiro, RJ, Brazil  Karin Wittmann W
<EMAIL_ADDRESS>Instituto de Física da UFRGS, Porto Alegre, RS, Brazil
Jorge Chávez-Carlos<EMAIL_ADDRESS>Department of Physics,
University of Connecticut, Storrs, Connecticut 06269, USA  Itzhak Roditi
<EMAIL_ADDRESS>Centro Brasileiro de Pesquisas Físicas/MCTI, 22290-180, Rio de
Janeiro, RJ, Brazil  Angela Foerster<EMAIL_ADDRESS>Instituto de Física
da UFRGS, Porto Alegre, RS, Brazil  Jorge G. Hirsch<EMAIL_ADDRESS>Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México,
Apdo. Postal 70-543, C.P. 04510 Cd. Mx., Mexico

###### Abstract

In this work, we investigate the semiclassical limit of a simple bosonic
quantum many-body system exhibiting both integrable and chaotic behavior. A
classical Hamiltonian is derived using coherent states. The transition from
regularity to chaos in classical dynamics is visualized through Poincaré
sections. Classical trajectories in phase space closely resemble the
projections of the Husimi functions of eigenstates with similar energy, even
in chaotic cases. It is demonstrated that this correlation is more evident
when projecting the eigenstates onto the Fock states. The analysis is carried
out at a critical energy where the eigenstates are maximally delocalized in
the Fock basis. Despite the imperfect delocalization, its influence is present
in the classical quantum properties under investigation. The study
systematically establishes quantum-classical correspondence for a bosonic
many-body system with more than two wells, even within the chaotic region.

##  I Introduction

As is well known, the field of quantum chaos begins by identifying singular
features present in quantum systems, obtained through the quantization of
chaotic classical systems [1, 2, 3, 4, 5]. While classical chaos is uniquely
associated with the explicit non-linear character of the equation of motion,
the term “quantum chaos” refers specifically to the system’s behavior in the
quantum regime. The most extensively studied example of this classical-quantum
correspondence is the billiard system, where precise and unambiguous relations
have been established [6, 7, 8]. When a chaotic system is quantized, the
separation between its eigenvalues exhibits a Wigner-Dyson distribution,
representing a typical signature for recognizing quantum chaos. Another
paradigmatic example with a well-established chaotic quantum-classical
correspondence is the Dicke model [9], and its study is a highly active
research field [10, 11, 12, 13, 14].

Other signatures are also employed to characterize standard features of
chaotic quantum systems, including correlation holes, a Gaussian distribution
of the density of states, and an off-diagonal distribution in the quantum
evolution of mean values of observables [15, 16, 17, 18, 19, 20, 21, 22].
These distinctive features are also present in quantum systems with no
explicit classical limit (such as nuclear or many-body quantum systems) and in
general random matrix models, where these signatures are understood as typical
universal behavior in interacting non-integrable quantum systems.

This work analyzes the quantum-classical correspondence of a three-well
bosonic model introduced in [23] and further studied in [24] from a pure
quantum chaos perspective. The model can be considered a generalized Bose-
Hubbard model with an integrable limit and exhibits mixed behavior in the
chaotic regime, as asserted by other recent studies of the Bose-Hubbard model
[25]. Perturbed integrable models with chaotic limits are intriguing
mathematical artifacts, as they provide the opportunity for a detailed study
of the integrable-chaotic transition, shedding light on the intricate aspects
of the mixing process. The studied model belongs to a family of integrable
bipartite models [26], 111For the four-well case, see [64] and [65], where the
integrable properties were studied to produce interferometry and NOON states.,
descending from particular two-well bosonic models. Integrability is
associated with nontrivial intrinsic symmetries determined by the Bethe
ansatz. Here, we thoroughly explore the quantum-classical correspondence in
this model for the integrable, mixed, and chaotic regimes.

The Bose-Hubbard model is intrinsically connected with the vast field of Bose-
Einstein condensates. It represents a specific limit of the spatial quantum
wave function within certain interacting potentials, describing the behavior
of ultracold atoms and diluted alkaline gases [27, 28]. Extensive literature
covers various aspects of this significant discipline, including localization,
catastrophe theory, quantum chaos, superfluidity, and the generation and
control of entanglement [29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55].

The paper is structured as follows: Section II defines the quantum
Hamiltonian, and we determine the classical critical energy
$E_{\mathrm{Classic}}$ through a pure eigenstates analysis. The nearest
eigenstates to $E_{\mathrm{Classic}}$ are also the more delocalized
eigenstates, and we select this energy region for our study. Subsequently, we
deduce the classical Hamiltonian written in an appropriate set of coordinates
and establish interesting properties in the classical trajectories. Section
III introduces coherent states and, from them, the Husimi function, which
measures quantum localization in phase space. The set of Fock states can also
be used for the same purpose. Remarkably, we find correspondence for many
classical trajectories in a short time and quantum localization in phase space
for many eigenstates in the chaotic regime. Specific patterns and
configurations are analyzed.

In Section IV, we study the previous classical-quantum correspondence, but now
for large-time classical trajectories. The study is conducted in three
regimes: regular, mixed, and chaotic. In the regular regime, we use the
associated symmetry to establish the correspondence in a precise form. In the
chaotic case, the correspondence is established using microcanonical mean
values of the Husimi functions, and we show that the classical pattern
obtained at large times is found in the quantum case because remanent
localization persists in the set of eigenstates. This is expressed as a slight
deviation from the Gaussian distribution of the components.

Section V contains the final discussion of the results obtained in this work.

##  II Quantum and classical Hamiltonians and the critical energy

###  II.1 Quantum Hamiltonian

The Hamiltonian of our model is an extension of a two-well potential version
studied in [56, 57]. It describes $N$ bosons in an aligned three-well
potential and is given by

$\displaystyle\hat{H}=$
$\displaystyle\frac{U}{N}\left(\hat{N}_{1}-\hat{N}_{2}+\hat{N}_{3}\right)^{2}+\epsilon\left(\hat{N}_{3}-\hat{N}_{1}\right)$
$\displaystyle+\frac{J}{\sqrt{2}}\left(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{2}^{\dagger}\hat{a}_{1}\right)+\frac{J}{\sqrt{2}}\left(\hat{a}_{2}^{\dagger}\hat{a}_{3}+\hat{a}_{3}^{\dagger}\hat{a}_{2}\right),$
(1)

where $\hat{N}_{k}=\hat{a}_{k}^{\dagger}\hat{a}_{k}$ is the number operator of
the $k$th-well, $\hat{a}_{k}^{\dagger}$ and $\hat{a}_{k}$ are respectively the
creation and annihilation operators of the $k$th-well. The total particle
number $N=N_{1}+N_{2}+N_{3}$ is fixed. $U$ is the magnitude of the boson
interactions, $J$ parameterize the jump between wells, and $\epsilon$ is an
external tilt. If $\epsilon=0$, the system is invariant by the interchange of
wells 1 and 3 and is integrable; if $\epsilon\neq 0$, the integrability is
broken.

The Hamiltonian (II.1) can be expressed in a matrix representation through
occupation number (Fock) states $|N_{1},N_{2},N_{3}\rangle$. The associated
Hamiltonian matrix has dimension $D\times D$ with $D=\frac{(N+2)!}{2!N!}$, and
its $D$ eigenvalues are obtained by numerical diagonalization.

The spectrum of this quantum Hamiltonian, for different parameter values,
exhibits both regularity and chaos regions. Various cases were analyzed in the
integrable limit $\epsilon=0$: the Rabi and Fock limits, which are highly
degenerate, and the Josephson regime in between. For $\epsilon>0$, a
significant number of avoided crossings are observed, and at $\epsilon\approx
1.5$, the spectrum displays the characteristic level repulsion of quantum
chaotic systems [24].

###  II.2 Classical critical energy and Shannon entropy

A comprehensive calculation of the equilibrium points associated with the
classical limit of the Hamiltonian (II.1) was conducted in [58]. The stable
configurations with minimal and maximal energy were identified, showing a
complete correspondence with the minimal and maximal eigenvalues for
sufficiently large $N$. While these results focus on ground-state properties,
the tools employed can also be utilized to compute intermediate energy
configurations, which, in principle, are associated with possible excited
quantum phase transitions.

These intermediate configurations are consistently present, particularly in
the vicinity of the chaotic regime identified by the parameters
$(U,J,\epsilon)=(0.7,1,1.5)$. In this configuration, an unstable critical
point arises with energy $E_{\mathrm{Classic}}\approx 0.0752$, represented by
the normalized state

$\mathbf{N}=(N_{1}/N,N_{3}/N,\phi_{12},\phi_{32})\approx(0.081,0.294,0,\pi),$

where $\phi_{12}$ and $\phi_{32}$ represent the phases in the classical
Hamiltonian approach, as detailed in the next subsection.

It has been found that the presence of classical critical points can be
reflected in the structure of quantum eigenstates. For instance, it was
demonstrated that the participation ratio is sensitive to the existence of the
critical point when considering only states in a basis with a particular
symmetry [59]. Regarding our model, the Shannon entropy can be examined for
each eigenstate in the Fock basis, and no evidence of abrupt changes
corresponding to the critical classical energy has been observed [24]. In Fig.
1(a), the Shannon entropy of the eigenstates is plotted against the
eigenenergy. The curved light blue line represents the mean value of the
Shannon entropy around the 200 eigenstates closest to each $E$. No singular
patterns are observed around the critical energy, marked with a vertical line.

We observe that the critical point with classical energy
$E_{\mathrm{Classic}}$ satisfies $N_{2}>N_{1}+N_{3}$ (From $\mathbf{N}$, we
see that $N_{2}/N\approx 0.625$). Considering only the Fock states with
$N_{2}>N_{1}+N_{3}$ (or equivalently with $N_{2}<N_{1}+N_{3}$), we can
determine the critical point energy using the Shannon entropy restricted to
this particular set. In Fig.1(b) it is shown an appreciable change in
$E/N\approx 0.0752$ when the Shannon entropy is restricted to the particular
set of Fock states that satisfy $N_{2}>N_{1}+N_{3}$ (orange points for
$N_{2}<N_{1}+N_{3}$).

Figure 1: We can express the total Shannon entropy as
$Sh_{\mathrm{Total}}^{E}=Sh_{N_{2}\geq
N_{1}+N_{3}}^{E}+Sh_{N_{2}<N_{1}+N_{3}}^{E}$; the total Shannon entropy is
depicted in purple in (a). In (b), the blue and orange figures correspond,
respectively, to $Sh_{N_{2}\geq N_{1}+N_{3}}^{E}$ and
$Sh_{N_{2}<N_{1}+N_{3}}^{E}$. The dashed blue line corresponds to the critical
energy, and the curved lines represent the mean value of the Shannon entropy
around the 200 eigenstates closest to each $E$. We use $\epsilon=1.5$ and
$N=120$.

As the total Shannon entropy reaches its maximum value at this point (see
Figure 1a), the components of the related eigenstates are more delocalized;
therefore, they are suitable for quantum-classical correspondence analysis in
the chaotic regime. For the analysis of the quantum chaotic regime, the
parameters $(U,J,\epsilon)=(0.7,1,1.5)$ will be adopted, and as a reference,
the energy $E_{\mathrm{Classic}}$. In quantum chaos, eigenstates with
eigenvalues $E/N$ closest to $E_{\mathrm{Classic}}$ will be considered. For
consistency, states of the mixed and regular (integrable) regimes will also be
analyzed using the same $(U,J)$ and energy parameters.

###  II.3 Classical Hamiltonian and dynamics

The classical Hamiltonian is obtained using coherent states,
$|\alpha\rangle=|\alpha_{1},\alpha_{2},\alpha_{3}\rangle$, where
${\alpha_{k}=\sqrt{N_{k}}\exp(i\phi_{k})}$. It leads to

$\displaystyle\mathcal{H}_{\text{cl}}$
$\displaystyle=\frac{\langle\alpha|\hat{H}|\alpha\rangle}{N}$ (2)
$\displaystyle=\frac{U}{N}\left(N_{1}-N_{2}+N_{3}\right)^{2}+\epsilon\left(N_{3}-N_{1}\right)$
$\displaystyle+$ $\displaystyle
J\sqrt{2}\left[\sqrt{N_{1}N_{2}}\cos(\phi_{1}-\phi_{2})+\sqrt{N_{2}N_{3}}\cos(\phi_{2}-\phi_{3})\right].$

Defining $\rho_{k}=\sqrt{N_{k}/N}$, the classical Hamiltonian reads

$\displaystyle\bar{\mathcal{H}}_{\text{cl}}=$
$\displaystyle\frac{\mathcal{H_{\text{cl}}}}{N}=U\left(\rho_{1}^{2}-\rho_{2}^{2}+\rho_{3}^{2}\right)^{2}+\epsilon\left(\rho_{3}^{2}-\rho_{1}^{2}\right)$
$\displaystyle+$ $\displaystyle
J\sqrt{2}\left[\rho_{1}\rho_{2}\cos(\phi_{1}-\phi_{2})+\rho_{2}\rho_{3}\cos(\phi_{2}-\phi_{3})\right]$
(3)
$\displaystyle=\frac{U}{4}\left(Q_{1}^{2}+P_{1}^{2}-Q_{2}^{2}-P_{2}^{2}+Q_{3}^{2}+P_{3}^{2}\right)^{2}$
$\displaystyle+\frac{\epsilon}{2}\left(Q_{3}^{2}+P_{3}^{2}-Q_{1}^{2}-P_{1}^{2}\right)$
$\displaystyle+\frac{J}{\sqrt{2}}\left[Q_{1}Q_{2}+P_{1}P_{2}+Q_{2}Q_{3}+P_{2}P_{3}\right]$
(4)

In the last line, we have introduced the canonical variables

$\displaystyle Q_{1}=\sqrt{2}\rho_{1}\cos(\phi_{1}),\,\,$ $\displaystyle
P_{1}=\sqrt{2}\rho_{1}\sin(\phi_{1}),$ (5) $\displaystyle
Q_{2}=\sqrt{2}\rho_{2}\cos(\phi_{2}),\,\,$ $\displaystyle
P_{2}=\sqrt{2}\rho_{2}\sin(\phi_{2}),$ (6) $\displaystyle
Q_{3}=\sqrt{2}\rho_{3}\cos(\phi_{3}),\,\,$ $\displaystyle
P_{3}=\sqrt{2}\rho_{3}\sin(\phi_{3})$ (7)

How many degrees of freedom has the above classical Hamiltonian? There are 3
generalized coordinates $Q_{1},Q_{2},Q_{3}$ and their conjugated momenta
$P_{1},P_{2},P_{3}$. There are two known conserved quantities: the energy $E$,
given that the Hamiltonian is time independent, and the total number of bosons
$N=N_{1}+N_{2}+N_{3}$, implying that
$\rho_{1}^{2}+\rho_{2}^{2}+\rho_{3}^{2}=1$.

It can also be seen in Eq. (II.1) and (II.3) that the Hamiltonian is also
invariant under the continuous transformation
$\phi_{k}\rightarrow\phi_{k}+\theta$, for an arbitrary phase $\theta$.

$\displaystyle Q^{\prime}_{k}=$
$\displaystyle\sqrt{2}\rho_{k}\cos(\phi_{k}+\theta)$ $\displaystyle=$
$\displaystyle\sqrt{2}\rho_{k}(\cos(\phi_{k})\cos(\theta)-\sin(\phi_{k})\sin(\theta)),$
$\displaystyle=$ $\displaystyle Q_{k}\cos(\theta)-P_{k}\sin(\theta),$ (8)
$\displaystyle\,\,P^{\prime}_{k}=$
$\displaystyle\sqrt{2}\rho_{k}\sin(\phi_{k}+\theta),$ $\displaystyle=$
$\displaystyle\sqrt{2}\rho_{k}(\sin(\phi_{k})\cos(\theta)+\cos(\phi_{k})\sin(\theta)),$
$\displaystyle=$ $\displaystyle P_{k}\cos(\theta)+Q_{k}\sin(\theta))$ (9)

It follows that
$Q_{k}^{2}-P_{k}^{2}={Q^{\prime}_{k}}^{2}-{P^{\prime}_{k}}^{2}$. Also

$\displaystyle Q^{\prime}_{1}Q^{\prime}_{2}+P^{\prime}_{1}P^{\prime}_{2}$
$\displaystyle=Q_{1}Q_{2}+P_{1}P_{2}$ (10)

Therefore, the Hamiltonian is invariant under this rotation in each phase
space by the same angle $\theta$. The trajectories $\\{Q_{k}(t),P_{k}(t)\\}$,
obtained as solutions of the dynamical Hamilton equations, will be described
by the same mathematical expressions in the different coordinate systems. It
implies that the trajectories with initial conditions
$\\{Q^{\prime}_{1}(\theta),Q^{\prime}_{2}(\theta),Q^{\prime}_{3}(\theta),P^{\prime}_{1}(\theta),P^{\prime}_{2}(\theta),P^{\prime}_{3}(\theta)\\}$
can be obtained from those with initial conditions
$\\{Q_{1}(0),Q_{2}(0),Q_{3}(0),P_{1}(0),P_{2}(0),P_{3}(0)\\}$ by applying the
transformation given by Eqs. (II.3).

With the above tools, we can consider expressing the Hamiltonian employing the
dynamical variables of the classical system, which are the $\rho_{k}$’s and
the phase differences $\phi_{k,k+1}=\phi_{k}-\phi_{k+1}$. Noticing that
$\rho_{2}=\sqrt{1-\rho_{1}^{2}-\rho_{3}^{2}}$ and that only the differences in
phases are relevant in the Hamiltonian, we can define
$\cos(\phi_{1}-\phi_{2})=\cos(\phi_{12})$ and
$\cos(\phi_{2}-\phi_{3})=\cos(\phi_{23})$.

The Hamiltonian takes the simplified form

$\displaystyle\bar{\mathcal{H}}_{\text{cl}}=$ $\displaystyle
U\left(2(\rho_{1}^{2}+\rho_{3}^{2})-1\right)^{2}+\epsilon\left(\rho_{3}^{2}-\rho_{1}^{2}\right)$
$\displaystyle+$ $\displaystyle
J\sqrt{2}\sqrt{1-\rho_{1}^{2}-\rho_{3}^{2}}\left[\rho_{1}\cos(\phi_{12})+\rho_{3}\cos(\phi_{23})\right]$
(11)

This Hamiltonian has 4 classical variables $\rho_{1},\phi_{12}$,
$\rho_{3},\phi_{23}$. We can introduce the canonical variables

$\displaystyle q_{1}=\sqrt{2}\rho_{1}\cos(\phi_{12}),\,\,$ $\displaystyle
p_{1}=\sqrt{2}\rho_{1}\sin(\phi_{12}),$ (12) $\displaystyle
q_{3}=\sqrt{2}\rho_{3}\cos(\phi_{23}),\,\,$ $\displaystyle
p_{3}=\sqrt{2}\rho_{3}\sin(\phi_{23}).$ (13)

In terms of these new variables, the Hamiltonian reads

$\displaystyle\bar{\mathcal{H}}_{\text{cl}}=U\left(q_{1}^{2}+p_{1}^{2}+q_{3}^{2}+p_{3}^{2}-1\right)^{2}$
$\displaystyle+\frac{\epsilon}{2}\left(q_{3}^{2}+p_{3}^{2}-q_{1}^{2}-p_{1}^{2}\right)$
$\displaystyle+J(q_{1}+q_{3})\sqrt{1-\frac{(q_{1}^{2}+p_{1}^{2}+q_{3}^{2}+p_{3}^{2})}{2}}.$
(14)

The classical Hamiltonian (4) has three classical degrees of freedom and two
conserved quantities: $E$ and $N$. Hamiltonian (II.3) has two degrees of
freedom and only one conserved quantity, the energy $E$. Both are equivalent
and display chaos, as we will show later.

The Hamilton equations have the form $\dot{x}_{i}=f_{i}(x_{1},\cdots,x_{m})$
for $i\in\\{1,\cdots m\\}$. In (II.3) the $f_{i}$ are simple polynomials in
the $x_{i}$’s variables and, in (II.3) the function $f_{i}$ are more involved.
It implies that in some cases, it is numerically more efficient to employ
(II.3) to calculate the classical trajectories using the Julia free software
Taylor.Integration [60].

The classical Hamiltonian depends on four variables. Throughout this work, we
choose $\rho_{2}=0$ as one of the initial conditions. This implies that
$\rho_{1}^{2}+\rho_{3}^{2}=1$, and the $J$ term in (II.3) is zero, indicating
that the classical energy is independent of the variables $\phi_{12}$ and
$\phi_{32}$. This condition never occurs when $\epsilon=0$, and for
$\epsilon>0$, the condition is possible for some energies. When the
$\rho_{2}=0$ point exists, all states with arbitrary phases $\phi_{12}$ and
$\phi_{32}$ have the same energy. In Figure 2, six trajectories are shown with
the initial condition $\rho_{2}=0$, fixing the same point in the plane
$(N_{1},N_{3})$, represented as a large blue dot. They have different phase
differences $\phi_{32}-\phi_{12}$, generating distinct trajectories with the
same energy.

Figure 2: Classical trajectories for $E_{\mathrm{Classic}}\approx 0.07524$ and
$(U,J,\epsilon)=(0.7,1,1.5)$ with an initial condition at $N_{2}=0$ (blue
point) and phase differences
$\phi_{32}-\phi_{12}\in{0,\pi/4,\pi/2,3\pi/2,\pi,5\pi/4}$, where
$\phi_{12}=-\pi$. The trajectories are distinct and originate from the same
$(N_{1},N_{3})$ value (the only one possible for $E_{\mathrm{Classic}}$).

Additionally, in the particular case where $\rho_{2}=0,\phi_{2}=0$, it is
possible to assume that $\phi_{12}=\phi_{1}$ and $\phi_{32}=\phi_{3}$. While
these two phases are rotated by the same angle $\theta$, this situation
differs from Eqs’s. (II.3) because the phase $\phi_{2}$ is held constant (and
is null). In this case, for different choices of $\phi_{12}$ and $\phi_{23}$
with equal $\phi_{32}-\phi_{12}$, the same trajectories are obtained but
rotated by an angle $\theta$ in the planes $(Q_{1},P_{1})$ and
$(Q_{3},P_{3})$. This effect can be observed in Figure 3, where three sets of
phases $\phi_{32},\phi_{12}$ are rotated by the same angle. In the plane
$(N_{1},N_{3})$, the trajectory remains the same, while in the planes
$(Q_{1},P_{1})$ and $(Q_{3},P_{3})$, the trajectories have the same shape but
are rotated.

Figure 3: Classical trajectories for $E_{\mathrm{Classic}}\approx 0.07524$ and
$(U,J,\epsilon)=(0.7,1,1.5)$ with the initial condition at $N_{2}=0$ (blue
point) and an equal phase difference $\phi_{32}-\phi_{12}=0$, with values
$\phi_{12}=\phi_{32}=0$, $\pi/4$, and $\pi/2$. The trajectories in the
$\left(N_{1},N_{3},\phi_{12},\phi_{32}\right)(t)$ coordinates are identical
(top figures a, b, and c), and the mutually equal changes in
$(\phi_{12},\phi_{32})$ rotate the trajectory in the $(Q_{1},P_{1})$ space
(middle figures) and the $(Q_{3},P_{3})$ space (bottom figures).

.

##  III Quantum localization in phase space

While classical trajectories cannot be directly mapped to quantum states due
to the uncertainty principle, the projection of a quantum state onto coherent
states allows us to obtain the Husimi quasi-probability distribution in phase
space.

The coherent state $|\alpha_{1},\alpha_{2},\alpha_{3}\rangle$ includes all
possible total numbers of bosons in the three wells. As the system under study
has a fixed total number of bosons, it is useful to introduce a coherent state
projected to a fixed total particle number, expressed as:

$\displaystyle|N_{1},$ $\displaystyle N_{3},\phi_{12},\phi_{32}\rangle=$
$\displaystyle\underset{n_{1}+n_{2}+n_{3}=N}{\sum_{n_{1},n_{2},n_{3}}}\left[P\left(n_{1},n_{2},n_{3},\frac{N_{1}}{N},\frac{N_{2}}{N},\frac{N_{3}}{N}\right)\right]^{1/2}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times
e^{in_{1}\phi_{12}}e^{in_{3}\phi_{32}}|n_{1},n_{2},n_{3}\rangle$ (15)

where $P$ is the multinomial distribution

$P\left(n_{1},n_{2},n_{3},p_{1},p_{2},p_{3}\right)=\frac{N!}{n_{1}!n_{2}!n_{3}!}p_{1}^{n_{1}}p_{2}^{n_{2}}p_{3}^{n_{3}},$
(16)

$\phi_{jk}=\phi_{j}-\phi_{k}$ and $p_{1}+p_{2}+p_{3}=1$. Note that to each
classical initial condition $(N_{1},N_{3},\phi_{12},\phi_{32})$ is associated
with the coherent state $|N_{1},N_{3},\phi_{12},\phi_{32}\rangle$.

The projection of the Husimi function of an eigenstate $|E\rangle$ with
eigenenergy $E$ over the plane $(N_{1},N_{3})$ can be calculated using

$\mathfrak{H}_{E}^{\mathcal{T}}(N_{1},N_{3})=\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\frac{d\phi_{12}}{2\pi}\frac{d\phi_{32}}{2\pi}\left|\langle
E|N_{1},N_{3},\phi_{12},\phi_{32}\rangle\right|^{2}$ (17)

An alternative description of the localization of the quantum states in the
same plane $(N_{1},N_{3})$ of the phase space is provided by the projection of
an eigenstate $|E\rangle$ on the Fock states

$\mathfrak{H}_{E}\left(N_{1},N_{3}\right)=|\langle
N_{1},N_{2}=N-N_{1}-N_{3},N_{3}|E\rangle|^{2}.$ (18)

The principal difference is that the quantity
$\mathfrak{H}_{E}\left(N_{1},N_{3}\right)$ generally exhibits a strong
interference pattern in non-chaotic regimes, with the advantage that the
computational numerical cost is very economical compared to the Husimi
function. Fock states have also been utilized in reference [61] to compare
quantum-classical correspondence, focusing on the quantum mean value evolution
in Fock states and caustics for the classical system at different times.

Figure 4: Projected Husimi function
$\mathfrak{H}_{E}^{\mathcal{T}}(N_{1},N_{3})$ (a) and the function
$\mathfrak{H}_{E}(N_{1},N_{3})$ using Fock states (b), for $E\approx
E_{\mathrm{Classic}}$, $\epsilon=0.7$ and $N=180$. We exponentiate each point
by 0.25 in both figures to obtain more details.

In Fig. 4, the two functions $\mathfrak{H}_{E}^{\mathcal{T}}(N_{1},N_{3})$ and
$\mathfrak{H}_{E}(N_{1},N_{3})$ are compared for an eigenstate with energy
close to $E_{\mathrm{Classic}}$, $N=150$, $\epsilon=0.7$, $U=0.7$ and $J=1$,
of the mixed regime. Note the similarity between the two figures, indicating
that in this particular projection, the function
$\mathfrak{H}_{E}(N_{1},N_{3})$ (which utilizes Fock states) provides similar
but sharper information than the Husimi function. In particular, the function
$\mathfrak{H}_{E}(N_{1},N_{3})$ in (b) exhibits a strong interference pattern,
offering better resolution than the Husimi function, with the additional
advantage of being very efficient in numerical computational terms compared to
the standard Husimi function.

###  III.1 Classical trajectories vs quantum projections on phase space

Here, we draw a parallel between classical trajectories and Husimi functions
for eigenvectors in the chaotic regime. In Figs 5(a,d,g,j,m,p), the
projections of six classical trajectories over the plane $(N_{1},N_{3})$ are
shown, evolved for a short time at energy $E_{\mathrm{Classic}}$ for
$\epsilon=1.5$, with six different initial conditions. They are compared with
the function $\mathfrak{H}_{E}(N_{1},N_{3})$ using Fock states in Figs
5(b,e,h,k,n,q), and the projected Husimi function
$\mathfrak{H}_{E}^{\mathcal{T}}(N_{1},N_{3})$ in Fig 5(c,f,i,l,o,r), for six
different eigenstates with eigenvalues near the classical energy
$E_{\mathrm{Classic}}$. Some more figures are in the appendix B.

Although the classical dynamics at energy $E_{\mathrm{Classic}}$ for
$\epsilon=1.5$ is chaotic, it is possible to observe particular trajectories
that are clearly recognized, where the density of points is larger, implying
that the system spends more time in these regions of the parameter space. They
could be associated with sticky orbits, like those found recently in chaotic
billiards [62]. A careful revision of the projections of the function
$\mathfrak{H}_{E}(N_{1},N_{3})$ and
$\mathfrak{H}_{E}^{\mathcal{T}}(N_{1},N_{3})$ belonging to hundreds of
eigenstates with eigenenergies close to $E_{\mathrm{Classic}}$ allowed to find
those whose projections closely resemble the classical trajectories. Similar
correspondence between the projections of the Husimi functions of eigenstates
have been reported for quasiperiodic orbits in the Dicke model of atoms in a
cavity [63, 11].

The quality of the projection over Fock states exceeds the one offered by the
Husimi functions. For this reason, only the first ones are employed in what
follows.

Figure 5: Classical trajectories projected onto the $(N_{1},N_{3})$ plane for
six distinct initial conditions with energy $E_{\mathrm{Classic}}$ (a, d, g,
j, m, p). Equivalent projections are obtained for the Husimi Functions for
some eigenstates with energy $E\approx E_{\mathrm{Classic}}$, projected into
both the Fock basis (18) (b, e, h, k, n, q) and onto coherent states (17) (c,
f, i, l, o, r). In all cases, $U=0.7$ and $\epsilon=1.5$. In quantum
projections, the number of particles varies: $N=130$ (b, k), $N=150$ (e, q),
and $N=270$ (h, n).

###  III.2 Trajectory pattern

The distinct shapes of the trajectories depicted in Fig. 5 were identified
through visual similarity analysis, spanning the entire energy spectrum of
systems with varying numbers of particles. These patterns emerge randomly for
different eigenstates, even within the chaotic quantum regime (II.1). Each
shape pattern can be found more than once within the chaotic regime of a
system, as shown in Fig. 6 for $N=210$. They arise even when the number of
particles in the system varies, as exemplified in Fig. 7, for $N=130$,
$N=180$, and $N=210$. It is also possible to identify the same pattern,
sparsely, in a wider range of the energy spectrum around the chaotic regime,
as exemplified in Fig. 8. Nonetheless, while it is possible to identify some
other sets of related patterns, only a few emerge reasonably well-defined.
Here, we present some of the most pronounced patterns (see Appendix A for some
more examples).

Figure 6: The function $\mathfrak{H}_{E}(N_{1},N_{3})$ (18) of three
eigenstates of a system with $N=210$, $U=0.7$, $\epsilon=1.5$ and $E\approx
E_{\mathrm{Classic}}$ that exhibit comparable projections. Figure 7: The
function $\mathfrak{H}_{E}(N_{1},N_{3})$ (18) for eigenstates of systems with
varying numbers of particles: (a) $N=130$, (b) $N=180$, and (c) $N=210$, with
energy $E\approx E_{\mathrm{Classic}}$. In all cases, $U=0.7$ and
$\epsilon=1.5$. Figure 8: The function $\mathfrak{H}_{E}(N_{1},N_{3})$ for six
specific eigenstates with eigenvalues (as indicated in the panels) distributed
over a wider range of the energy spectrum, around $E_{\mathrm{class}}$. The
distributions tend to be skewed, but it is possible to relate them to a
standard shape. In all cases, $N=270$, $U=0.7$ and $\epsilon=1.5$.

The intriguing observation is that these patterns appear recurrently, also for
classical trajectories of different initial conditions within the chaotic
regime. This can be verified, for example, in Fig. 2. In this context, these
projection patterns reflect a characteristic that pervades both classical and
quantum systems within the bosonic model under investigation. This subject
will be explored further in the next section.

##  IV The classical-quantum correspondence

This section explores the quantum-classical correspondence for the 3-well
model in the regular ($\epsilon=0$), mixed ($\epsilon=0.7$), and chaotic
($\epsilon=1.5$) regimes, with $U=0.7$ and $J=1$. The analysis will center
around the energy $E_{\mathrm{Classic}}\approx 0.0752$ (the critical point
energy studied in ref. [58]), as in the preceding sections.

For the classical dynamics, the plots $N_{1}/N$ vs $\phi_{12}$ on the plane
$\phi_{32}=0$ are Poincaré sections, useful to visualize the regular, mixed or
chaotic trajectories. These sections are constructed in two ways: using
hundreds of initial conditions evolved for short times and a few initial
conditions evolved for long times. They provide similar yet complementary
views of classical dynamics. Some of these initial conditions are evolved in
time to obtain projected trajectories in coordinates $N_{3}/N$ vs $N_{1}/N$,
which are then compared with the quantum probability distributions
$\mathfrak{H}_{E}(N_{1},N_{3})$.

###  IV.1 The integrable case

We select the parameter values $\epsilon=0$, $U=0.7$, and $J=1$ to study the
integrable case.

####  IV.1.1 The conserved quantity $\hat{Q}$

The quantum operator $\hat{Q}$

$\hat{Q}=\left(\hat{N}_{1}+\hat{N}_{3}\right)-\left(a_{1}^{\dagger}a_{3}+a_{3}^{\dagger}a_{1}\right)$
(19)

is in this case a constant of motion, satisfying
$\left[\hat{Q},\hat{H}\right]=0$. $\hat{Q}$ can be obtained directly using the
Bethe ansatz [23]. There are $N+1$ eigenvalues of $\hat{Q}/N$, equally
distributed in the interval $[0,2]$ with $q_{1}=0$, $q_{N+1}=2$ and
$q_{m}=2(m-1)/N$. Each eigenstate $|E\rangle$ is also an eigenvector of
$\hat{Q}/N$. Different eigenstates have the same eigenvalues $q$ of
$\hat{Q}/N$, with a degeneracy $m$-dependent: $N+2-m$.

The classical limit of $\hat{Q}$ in normalized coordinates
$\rho_{k}=\sqrt{N_{k}/N}$ is

$\frac{Q}{N}=\rho_{1}^{2}+\rho_{3}^{2}-2\rho_{1}\rho_{3}\cos\left(\phi_{32}-\phi_{12}\right)$
(20)

$Q/N$ takes all the possible values in the interval $\left[0,2\right]$.
Particularly, $Q/N=0$ when $N_{2}/N=1$ and $Q/N=2$ when $N_{1}/N=N_{3}/N=1/2$,
$\phi_{32}-\phi_{12}=\pm n\pi$, with $n\in\mathbb{Z}$.

Two classical conserved quantities exist: the total energy and the classical
$Q$, implying two natural frequencies. This work focuses on quantum-classical
correspondence, so we do not calculate the frequencies. However, the evolution
of many initial conditions suggests that the frequencies are incommensurable
and manifest in a quasi-periodic motion.

The quantum-classical correspondence is straightforward in this case: any
eigenstate with eigenvalue $q$ of $\hat{Q}/N$ has a function
$\mathfrak{H}_{E}(N_{1},N_{3})$ that is essentially identical to a projected
classical trajectory with the classical value $Q/N=q$ and a classical energy
value very close to $E$.

####  IV.1.2 Classical trajectories

The integrability and quasiperiodicity allow for the classification of all
classical trajectories. For each classical value of $Q/N$, there is a unique
trajectory encompassing a specific region. Figure 9 b) displays the
intersection of six trajectories with six different values of $Q$ and the
hypersurface $\phi_{32}=0$. The points are projected onto the coordinates
$(N_{1},\phi_{12})$. Well-localized bands are observed, associated with
quasiperiodic orbits.

Fig.9 a) is a similar Poincaré section, for many trajectories with initial
conditions distributed in all the hypersurface evolved for short times. The
patterns are similar, but the bands are harder to observe.

Figure 9: a) Coordinates $(N_{1},\phi_{12})$ of the Poincaré section
$\phi_{32}=0$, constituted by 460 different initial conditions with energy
$E=E_{\mathrm{Classic}}$. b)$(N_{1},\phi_{12})$ coordinates of the
intersection points of six trajectories (with different values of $Q$ equal to
$0.1,0.4,0.8,1.0,1.3,1.5$ respectively and energy $E_{\mathrm{Classic}}$) and
the 2D surface $(E,\phi_{32})=(E_{\mathrm{Classic}},0)$ evolved in large times
($t_{f}=100000$).

To explore the quantum-classical correspondence, the classical trajectories at
long times projected in the coordinates $(N_{1},N_{3})$ are compared with the
quantum projection $\mathfrak{H}_{E}(N_{1},N_{3})$.

Figure 10:  Upper panels (a, b, c): Function $\mathfrak{H}_{E}(N_{1},N_{3})$
for three eigenstates with eigenvalues closest to $E_{\mathrm{Classic}}$ and
different values of $q$ (0.1, 0.4, and 1, respectively) of $\hat{Q}$. Bottom
panels (d, e, f): Projected classical trajectories on the plane
$(N_{1},N_{3})$ for classical initial conditions with $Q$ values of $0.1$,
$0.4$, and $1.0$.

Figure 10 is an example of the quantum-classical correspondence for
$\epsilon=0$. The function $\mathfrak{H}_{E}(N_{1},N_{3})$ for three
eigenstates with energies closest to $E_{\mathrm{Classic}}$ and three
different values of $q$ (respectively 0.1, 0.4, and 1, from left to right) is
compared with three classical trajectories with the same value of $Q/N$ and
classical energy $E_{\mathrm{Classic}}$ projected onto the plane
$(N_{1},N_{3})$. The remarkable similarity between both descriptions and the
fringes resembling an interference pattern are also visible in the classical
plots.

###  IV.2 The mixed case $\epsilon=0.7$

In this section, we study the case $\epsilon=0.7$, which shows a mixed
behavior.

The coexistence of regular and chaotic regions can be observed in the Poincaré
sections shown in Fig. 11, displaying the results for the fixed energy
$E_{\mathrm{Classic}}$ and $\phi_{32}=0$. Fig. 11 (a) was made using 513
initial conditions evolved at short times, and Fig. 11 (b) for five initial
conditions at long times. The regions of regular dynamics with quasi-periodic
trajectories and the chaotic ones are recognized in both figures. The purple
points set up a typical chaotic region, and the other four trajectories
configure well-localized tiny strips.

Figure 11: Poincaré sections for $\epsilon=0.7$, $E=E_{\mathrm{Classic}}$ and
$\phi_{32}=0$. a) $513$ trajectories homogeneously distributed and evolved for
short times. b) Five initial conditions evolved for long times.

For $\epsilon=0.7$, analyzing the quantum-classical correspondence is more
challenging. On the classical side, the five trajectories in Fig. 11 (for all
times) are projected onto the coordinates $(N_{1},N_{3})$. Three are presented
at the bottom of Fig. 12. The trajectories are projected onto the coordinates
$(N_{1},N_{3})$ for times from $t=0$ to $t=10000$. The purple classical
trajectory is chaotic, and the other two are regular and quasiperiodic. For
the quantum analysis, 200 eigenstates closest to $E_{\mathrm{Classic}}$
calculated with $N=180$ were selected, and their functions
$\mathfrak{H}_{E}(N_{1},N_{3})$ were calculated. Most of the 200 eigenstates
extend in the same region as the purple trajectory; only one is presented. A
few eigenstates in this energy region have a localized function
$\mathfrak{H}_{E}(N_{1},N_{3})$, resembling the classical projections. Two of
them are shown in Fig. 12(b) and 12(c).

Figure 12: Upper panels (a, b, c): Function $\mathfrak{H}_{E}(N_{1},N_{3})$
for three eigenstates with eigenvalues closest to $E_{\mathrm{Classic}}$.
Bottom panels (d, e, f): Projected classical trajectories onto the plane
$(N_{1},N_{3})$ for initial conditions with $E=E_{\mathrm{Classic}}$ and
$\epsilon=0.7$.

###  IV.3 The chaotic case $\epsilon=1.5$

The study of the case with $(U,J,\epsilon)=(0.7,1,1.5)$ reveals explicit
quantum chaos behavior in the middle of the energy spectrum region
(particularly at $E/N\approx 0.0752$), and in this case, the associated
classical scenario exhibits a typical chaotic Poincaré section (see Fig. 13).
For any initial condition with energy $E=E_{\mathrm{Classic}}\approx 0.0752$,
each trajectory encompasses all accessible phase space and has the same long-
time structure.

Figure 13: Poincaré sections in coordinates $(N_{1}/N,N_{3}/N)$ for
$(E,\phi_{32})=(E_{\mathrm{Classic}},\pi)$ (a) and
$(E,\phi_{32})=(E_{\mathrm{Classic}},0)$ (b). An overlap of both (a) and (b)
Poincaré sections is shown in (c). Panel (d) displays a Poincaré section for
an arbitrary classical chaotic trajectory for long times. The ’x’ points mark
the classical unstable critical point, while the magenta circles mark the
prevailing initial state of the classical trajectories.

Fig. 13(a) displays the Poincaré section on $(N_{1}/N,N_{3}/N)$ with energy
$E_{\mathrm{Classic}}$, but restricted to $\phi_{32}=\pi$, obtained from one
initial condition at long times. It covers two areas connected at one point,
which coincides with the unstable equilibrium points studied in [58]. Fig.
13(b) displays the Poincaré section restricted to $\phi_{32}=0$. Both exhibit
a noticeable concentration of points in their external contour, where the
values of $\phi_{12}$ are also very close to ${0,\pm\pi}$.

As shown in ref. [58], any classical critical point has exactly these extremal
phase values. As a consequence, any chaotic trajectory spends more time close
to these contour points since, in them,
$\partial\mathcal{H}/\partial\phi_{i2}=0$ with $i\in{1,3}$. Fig. 13d) displays
the $(N_{1},N_{3})$ projection of an evolved trajectory from an arbitrary
initial condition. The density of points is higher in the vicinity of the
critical regions.

To analyze the quantum description of the phase space in the chaotic regime,
it is useful to calculate the average of the distributions
$\mathfrak{H}_{E}(N_{1},N_{3})$ in an energy window $\mathfrak{E}$ in the
interval $\left(E_{\mathrm{Classic}}-\Delta E/2,E_{\mathrm{Classic}}+\Delta
E/2\right)$:

$\bar{\mathfrak{H}}(N_{1},N_{3})=\frac{1}{\mathcal{N}}\sum_{E\in\mathfrak{E}}\mathfrak{H}_{E}(N_{1},N_{3})$
(21)

It is shown for $N=180$, with $\mathfrak{E}$ containing 172 eigenstates around
$E/N\approx 0.07523767$ with $\Delta E=0.02$ in Fig. 14.

Figure 14: Panel (a) displays an arbitrary classical chaotic trajectory with
an energy of $E/N\approx 0.07523767$. In panel (b), the mean value
distribution, represented by Eq. (21), is presented. This distribution
encompasses all eigenstates falling within an energy window of $\Delta E=0.02$
(comprising 172 eigenstates) centered at $E_{\mathrm{Classic}}\approx
0.07523767$ for a system of size $N=180$. It is noteworthy to observe the
remarkable correspondence between the density distributions of the classical
and quantum projections.

It is remarkable that the average over a small energy window of the projection
of the quantum eigenstates closely reproduces the phase space structure
observed for any chaotic trajectory when projected onto the coordinates
$(N_{1},N_{3})$. In this chaotic regime, the long-time behavior of any
classical trajectory with classical energy $E_{\mathrm{Classic}}$ is reflected
in the microcanonical average of the projected Husimi function around the
nearest $\mathcal{N}$’s eigenenergies to $E_{\mathrm{Classic}}$. Any classical
chaotic trajectory spends more time in the dense region shown in Fig. 14,
which is directly connected with the unstable critical point with energy
$E/N\approx 0.0752$ having the same values of
$(\phi_{12},\phi_{32})=(0,\pm\pi)$.

Figure 15: The mean value distribution, given by Eq. (21), is showcased across
various energy windows, each centered at $E_{\mathrm{Classic}}\approx
0.07523767$. These windows correspond to different numbers of eigenstates: (a)
1, (b) 10, (c) 50, (d) 100, (e) 200, and (f) 300, for a system of size
$N=180$.

In Fig. 15, we display the mean value (21) for different values of $\Delta E$
with $N=180$ and eigenvalues $E/N$ closest to $E_{\mathrm{Classic}}$. Each
energy interval contains a different number of associated eigenvectors. It is
clear from the figure that the characteristic pattern in the average (21) is
generated for different $\mathcal{N}$, even present when $\mathfrak{E}$
contains only ten eigenvectors (Fig. 15b). Fig. 15e shows a more fine
contrast, having $\mathcal{N}=200$. On the other hand, for $\mathcal{N}=300$,
the pattern begins to defocus, mainly in the region around the critical point.

Maximally chaotic eigenstates exhibit a specific Gaussian distribution of
their components; therefore, mean values like (21) would eliminate any
specific pattern of the particular eigenstate distributions. The more chaotic
eigenstates do not follow a perfectly Gaussian distribution of the components,
although the system satisfies all the spectrum criteria of chaotic quantum
systems [24]. Particularly, the deviation from Gaussian behavior is more
evident at the extremum of the distribution, where the more prominent
components reside. This can be interpreted as a remnant localization, which is
not apparent for individual eigenstates but is evident in the mean value (21)
where many close eigenvalues are present. The dense region in figure 14b
expresses the slight deviation from the total delocalization characteristic of
chaotic eigenstates. The Fock states in the dense region contribute more to
the eigenstate superposition. This is visible in Fig. 16, where we show the 20
Fock states with the largest absolute value contribution in each one of the
$\mathcal{N}=200$ eigenstates used in (21). The 4000 Fock states, with the
principal contribution, constitute the dense region of Fig. 14b.

Figure 16:  $(N_{1}/N,N_{3}/N)$ indices of the Fock states with the most
prominent components in the 200 eigenstates used in Eq. (21). The total number
of particles in each case is (a) $N=100$, (b) $N=250$, and (c) $N=500$. From
each eigenstate, 20 components were selected. It can be seen that the
classical and quantum dense set region in Fig. 14 is reproduced.

##  V Discussion

In this work, we have established the quantum-classical correspondence of a
many-body bosonic system with an integrable limit. The model is sufficiently
simple to apply many techniques of the quantum-chaos discipline but has some
particularities, such as the quasiperiodic behavior of the classical
trajectories and the non-standard distribution of the chaotic states. The
transition to chaos is demonstrated, and a detailed comparison between
classical and quantum is carried out. In the integrable limit, a complete
classification of the orbits is recognized, and the associated projected
Husimi function of an eigenstate with energy very close to the classical
energy is calculated. In particular, we find an eigenstate with energy very
close to $E_{\mathrm{Classic}}$ for sufficiently large $N$ associated with any
classical quasiperiodic trajectory. The transition to chaos mixes the
quasiperiodic orbits and the associated quantum states in a complicated way;
the quantum states are now related to intricate trajectories only in short
time intervals. Remarkably, the large-time structure of any chaotic trajectory
is reproduced by the averaged Husimi function, which represents a dense region
connected with the unstable critical point that characterizes this larger time
structure. The dense region is also related to the Fock states with the larger
contribution in each one of the eigenstates belonging to $\mathfrak{E}$, used
in the Husimi mean value (21). Our study demonstrates how a slight deviation
from perfectly delocalized eigenstates influences the large-time behavior of
classical chaotic trajectories and establishes the quantum-classical
correspondence firmly, even in the chaotic case.

##  VI acknowledgment

E. Castro thanks the Brazilian CNPq agency for partial financial support. AF
and KW acknowledge support from CNPq (Conselho Nacional de Desenvolvimento
Científico e Tecnológico) - Edital Universal 406563/2021-7. J.G.H acknowledge
partial financial support from project PAPIIT-UNAM IN109523. We use the Julia
free software TaylorIntegration.jl for the dynamics of all the classical
trajectories. We thank Lea Ferreira dos Santos for valuable discussions and
suggestions.

##  Appendix A Quantum projection patterns: Additional figures

Here, we present additional examples of distinct forms that prominently emerge
in the projections of the functions,
$\mathfrak{H}_{E}\left(N_{1},N_{3}\right)$, for quantum chaotic systems. These
projections correspond to eigenvectors with energies close to the value of the
unstable classical critical point, $E\approx E_{\mathrm{Classic}}$. In Fig.
17, two recurring shape patterns are highlighted, emerging through different
eigenvectors from the most chaotic region of the quantum system spectrum, with
a fixed $N=270$. Conversely, in Fig. 18, shape patterns resulting from
eigenvectors of systems with different numbers of particles are presented.

Figure 17: The function $\mathfrak{H}_{E}(N_{1},N_{3})$ of two particular sets
of eigenstates that exhibit similar distributions, with $E\approx
E_{\mathrm{Classic}}$, for $N=270$, $U=0.7$ and $\epsilon=1.5$. Figure 18: The
function $\mathfrak{H}_{E}(N_{1},N_{3})$ of a set of eigenstates, with
$E\approx E_{\mathrm{Classic}}$, that display similar distributions for
systems with different numbers of particles, respectively: (upper) $N=180$,
$N=220$ and $N=270$. (middle) $N=130$, $N=160$ and $N=260$. (bottom) $N=180$,
$N=200$ and $N=240$. For all, $U=0.7$ and $\epsilon=1.5$.

##  Appendix B Classical-quantum correspondence patterns: Additional figures

It is possible to identify classical-quantum correspondence for all the
standard forms shown in this article. To strengthen the argument, in Fig. 19
are shown three additional examples, showcasing projection forms observed in
both the classical trajectories of with critical energy $E_{\mathrm{Classic}}$
and in the quantum Husimi functions
$\mathfrak{H}_{E}^{\mathcal{T}}(N_{1},N_{3})$ and
$\mathfrak{H}_{E}(N_{1},N_{3})$, for eigenstates with eigenvalues $E\approx
E_{\mathrm{Classic}}$. Identification is facilitated through visual
comparison.

Figure 19: Projections of three classical trajectories onto the plane
$(N_{1},N_{3})$ for different initial conditions with energy
$E_{\mathrm{Classic}}$ (a, d, g). Equivalent projections are obtained for the
function $\mathfrak{H}_{E}(N_{1},N_{3})$ (b, e, h) and
$\mathfrak{H}_{E}^{\mathcal{T}}(N_{1},N_{3})$ (c, f , i), for eigenstates with
energy $E\approx E_{\mathrm{Classic}}$. In quantum projections, the number of
particles varies: $N=130$ (top row), $N=200$ (middle row), and $N=270$ (bottom
row). In all cases, $U=0.7$, $\epsilon=1.5$, and $J=1$.

## References

  * [1] F. Haake, Quantum Signatures of Chaos.  Berlin, Heidelberg: Springer-Verlag, 2010.
  * [2] M. C. Gutzwiller, Chaos in Classical and Quantum Mechanics.  New York, NY: Springer, 2013.
  * [3] M. C. Gutzwiller, “Periodic Orbits and Classical Quantization Conditions,” Journal of Mathematical Physics, vol. 12, pp. 343–358, 10 2003.
  * [4] R. Aurich, J. Bolte, and F. Steiner, “Universal signatures of quantum chaos,” Phys. Rev. Lett., vol. 73, pp. 1356–1359, Sep 1994.
  * [5] O. Bohigas, S. Tomsovic, and D. Ullmo, “Manifestations of classical phase space structures in quantum mechanics,” Physics Reports, vol. 223, no. 2, pp. 43–133, 1993.
  * [6] O. Bohigas, M. J. Giannoni, and C. Schmit, “Characterization of chaotic quantum spectra and universality of level fluctuation laws,” Phys. Rev. Lett., vol. 52, pp. 1–4, Jan 1984.
  * [7] P. Sěba, “Wave chaos in singular quantum billiard,” Phys. Rev. Lett., vol. 64, pp. 1855–1858, Apr 1990.
  * [8] C. Lozej, G. Casati, and T. Prosen, “Quantum chaos in triangular billiards,” Phys. Rev. Res., vol. 4, p. 013138, Feb 2022.
  * [9] R. H. Dicke, “Coherence in spontaneous radiation processes,” Phys. Rev., vol. 93, pp. 99–110, Jan 1954.
  * [10] L. Bakemeier, A. Alvermann, and H. Fehske, “Dynamics of the dicke model close to the classical limit,” Phys. Rev. A, vol. 88, p. 043835, Oct 2013.
  * [11] S. Pilatowsky-Cameo, D. Villaseñor, M. A. Bastarrachea-Magnani, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, “Quantum scarring in a spin-boson system: fundamental families of periodic orbits,” New Journal of Physics, vol. 23, p. 033045, mar 2021.
  * [12] C. Emary and T. Brandes, “Chaos and the quantum phase transition in the dicke model,” Phys. Rev. E, vol. 67, p. 066203, Jun 2003.
  * [13] M. A. Bastarrachea-Magnani, B. L. del Carpio, S. Lerma-Hernández, and J. G. Hirsch, “Chaos in the dicke model: quantum and semiclassical analysis,” Physica Scripta, vol. 90, p. 068015, may 2015.
  * [14] Q. Wang, “Quantum chaos in the extended dicke model,” Entropy, vol. 24, no. 10, 2022.
  * [15] T. Gorin and T. H. Seligman, “Signatures of the correlation hole in total and partial cross sections,” Phys. Rev. E, vol. 65, p. 026214, Jan 2002.
  * [16] E. J. Torres-Herrera and L. F. Santos, “Dynamical manifestations of quantum chaos: correlation hole and bulge,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 375, no. 2108, p. 20160434, 2017.
  * [17] L. F. Santos, F. Borgonovi, and F. M. Izrailev, “Onset of chaos and relaxation in isolated systems of interacting spins: Energy shell approach,” Phys. Rev. E, vol. 85, p. 036209, Mar 2012.
  * [18] J. French and S. Wong, “Validity of random matrix theories for many-particle systems,” Physics Letters B, vol. 33, no. 7, pp. 449–452, 1970.
  * [19] 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., vol. 53, pp. 385–479, Jul 1981.
  * [20] L. F. Santos and M. Rigol, “Onset of quantum chaos in one-dimensional bosonic and fermionic systems and its relation to thermalization,” Phys. Rev. E, vol. 81, p. 036206, Mar 2010.
  * [21] R. Mondaini and M. Rigol, “Eigenstate thermalization in the two-dimensional transverse field ising model. ii. off-diagonal matrix elements of observables,” Phys. Rev. E, vol. 96, p. 012157, Jul 2017.
  * [22] L. F. Santos, F. Pérez-Bernal, and E. J. Torres-Herrera, “Speck of chaos,” Phys. Rev. Res., vol. 2, p. 043034, Oct 2020.
  * [23] K. W. Wilsmann, L. H. Ymai, A. P. Tonel, J. Links, and A. Foerster, “Control of tunneling in an atomtronic switching device,” Commun Phys, vol. 1, p. 91, Dec 2018.
  * [24] K. Wittmann W., E. R. Castro, A. Foerster, and L. F. Santos, “Interacting bosons in a triple well: Preface of many-body quantum chaos,” Phys. Rev. E, vol. 105, p. 034204, Mar 2022.
  * [25] G. Nakerst and M. Haque, “Chaos in the three-site bose-hubbard model: Classical versus quantum,” Phys. Rev. E, vol. 107, p. 024210, Feb 2023\.
  * [26] L. H. Ymai, A. P. Tonel, A. Foerster, and J. Links, “Quantum integrable multi-well tunneling models,” J. Phys. A Math. Theor., vol. 50, no. 26, p. 264001, 2017.
  * [27] A. J. Leggett, “Bose-einstein condensation in the alkali gases: Some fundamental concepts,” Rev. Mod. Phys., vol. 73, pp. 307–356, Apr 2001\.
  * [28] I. Bloch, J. Dalibard, and W. Zwerger, “Many-body physics with ultracold gases,” Rev. Mod. Phys., vol. 80, pp. 885–964, Jul 2008.
  * [29] T. F. Viscondi and K. Furuya, “Dynamics of a bose–einstein condensate in a symmetric triple-well trap,” Journal of Physics A: Mathematical and Theoretical, vol. 44, p. 175301, mar 2011.
  * [30] L. Cao, I. Brouzos, S. Zöllner, and P. Schmelcher, “Interaction-driven interband tunneling of bosons in the triple well,” New Journal of Physics, vol. 13, p. 033032, mar 2011.
  * [31] G. M. Koutentakis, S. I. Mistakidis, and P. Schmelcher, “Quench-induced resonant tunneling mechanisms of bosons in an optical lattice with harmonic confinement,” Phys. Rev. A, vol. 95, p. 013617, Jan 2017.
  * [32] I. Bloch, “Ultracold quantum gases in optical lattices,” Nat Phys, vol. 1, pp. 23–30, Oct 2005.
  * [33] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, “Boson localization and the superfluid-insulator transition,” Phys. Rev. B, vol. 40, pp. 546–570, Jul 1989.
  * [34] C. Kollath, G. Roux, G. Biroli, and A. M. Läuchli, “Statistical properties of the spectrum of the extended bose–hubbard model,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2010, p. P08011, aug 2010\.
  * [35] S. Dutta, M. C. Tsatsos, S. Basu, and A. U. J. Lode, “Management of the correlations of ultracoldbosons in triple wells,” New Journal of Physics, vol. 21, p. 053044, may 2019.
  * [36] N. Oelkers and J. Links, “Ground-state properties of the attractive one-dimensional bose-hubbard model,” Phys. Rev. B, vol. 75, p. 115119, Mar 2007.
  * [37] Q. Guo, X. Chen, and B. Wu, “Tunneling dynamics and band structures of three weakly coupled bose-einstein condensates,” Opt. Express, vol. 22, pp. 19219–19234, Aug 2014.
  * [38] A. R. Kolovsky and A. Buchleitner, “Quantum chaos in the bose-hubbard model,” Europhysics Letters, vol. 68, p. 632, nov 2004.
  * [39] T. Choy and F. Haldane, “Failure of bethe-ansatz solutions of generalisations of the hubbard chain to arbitrary permutation symmetry,” Physics Letters A, vol. 90, no. 1, pp. 83–84, 1982.
  * [40] A. I. Streltsov, K. Sakmann, O. E. Alon, and L. S. Cederbaum, “Accurate multi-boson long-time dynamics in triple-well periodic traps,” Phys. Rev. A, vol. 83, p. 043604, Apr 2011.
  * [41] T. Lahaye, T. Pfau, and L. Santos, “Mesoscopic ensembles of polar bosons in triple-well potentials,” Phys. Rev. Lett., vol. 104, p. 170404, Apr 2010\.
  * [42] P. Buonsante, R. Franzosi, and V. Penna, “Control of unstable macroscopic oscillations in the dynamics of three coupled bose condensates,” Journal of Physics A: Mathematical and Theoretical, vol. 42, p. 285307, jun 2009\.
  * [43] A. Foerster and E. Ragoucy, “Exactly solvable models in atomic and molecular physics,” Nuclear Physics B, vol. 777, no. 3, pp. 373–403, 2007.
  * [44] A. P. Tonel, L. H. Ymai, K. W. W., A. Foerster, and J. Links, “Entangled states of dipolar bosons generated in a triple-well potential,” SciPost Phys. Core, vol. 2, p. 3, 2020.
  * [45] G. Nakerst and M. Haque, “Eigenstate thermalization scaling in approaching the classical limit,” Phys. Rev. E, vol. 103, p. 042109, Apr 2021.
  * [46] M. Rautenberg and M. Gärttner, “Classical and quantum chaos in a three-mode bosonic system,” Phys. Rev. A, vol. 101, p. 053604, May 2020.
  * [47] S. Ray, D. Cohen, and A. Vardi, “Chaos-induced breakdown of bose-hubbard modeling,” Phys. Rev. A, vol. 101, p. 013624, Jan 2020.
  * [48] S. Bera, R. Roy, A. Gammal, B. Chakrabarti, and B. Chatterjee, “Probing relaxation dynamics of a few strongly correlated bosons in a 1d triple well optical lattice,” Journal of Physics B: Atomic, Molecular and Optical Physics, vol. 52, p. 215303, oct 2019.
  * [49] A. Richaud and V. Penna, “Phase separation can be stronger than chaos,” New Journal of Physics, vol. 20, p. 105008, oct 2018.
  * [50] M. A. Garcia-March, S. van Frank, M. Bonneau, J. Schmiedmayer, M. Lewenstein, and L. F. Santos, “Relaxation, chaos, and thermalization in a three-mode model of a bose–einstein condensate,” New Journal of Physics, vol. 20, p. 113039, nov 2018.
  * [51] L. Guo, L. Du, C. Yin, Y. Zhang, and S. Chen, “Dynamical evolutions in non-hermitian triple-well systems with a complex potential,” Phys. Rev. A, vol. 97, p. 032109, Mar 2018.
  * [52] G. McCormack, R. Nath, and W. Li, “Hyperchaos in a bose-hubbard chain with rydberg-dressed interactions,” Photonics, vol. 8, no. 12, 2021.
  * [53] T. Yan, M. Collins, R. Nath, and W. Li, “Signatures of quantum chaos of rydberg-dressed bosons in a triple-well potential,” Atoms, vol. 11, no. 6, 2023.
  * [54] C. B. Dağ, S. I. Mistakidis, A. Chan, and H. R. Sadeghpour, “Many-body quantum chaos in stroboscopically-driven cold atoms,” Commun Phys, vol. 6, p. 136, Jun 2023.
  * [55] K. W. W., L. H. Ymai, B. H. C. Barros, J. Links, and A. Foerster, “Controlling entanglement in a triple-well system of dipolar atoms,” Phys. Rev. A, vol. 108, p. 033313, Sep 2023.
  * [56] A. P. Tonel, J. Links, and A. Foerster, “Quantum dynamics of a model for two Josephson-coupled Bose–Einstein condensates,” Journal of Physics A: Mathematical and General, vol. 38, p. 1235, jan 2005.
  * [57] J. Links, A. Foerster, A. Tonel, and G. Santos, “The two-site Bose-Hubbard model,” Ann. Henri Poincaré, vol. 7, p. 1591, 2006.
  * [58] E. R. Castro, J. Chávez-Carlos, I. Roditi, L. F. Santos, and J. G. Hirsch, “Quantum-classical correspondence of a system of interacting bosons in a triple-well potential,” Quantum, vol. 5, p. 563, Oct. 2021.
  * [59] L. F. Santos, M. Távora, and F. Pérez-Bernal, “Excited-state quantum phase transitions in many-body systems with infinite-range interaction: Localization, dynamics, and bifurcation,” Phys. Rev. A, vol. 94, p. 012113, Jul 2016.
  * [60] J. A. Pérez-Hernández and L. Benet, “Perezhz/taylorintegration.jl:,” 10.5281/zenodo.2562352, vol. v0.4.1,, 2019.
  * [61] W. Kirkby, Y. Yee, K. Shi, and D. H. J. O’Dell, “Caustics in quantum many-body dynamics,” Phys. Rev. Res., vol. 4, p. 013105, Feb 2022.
  * [62] i. c. v. Lozej, D. Lukman, and M. Robnik, “Effects of stickiness in the classical and quantum ergodic lemon billiard,” Phys. Rev. E, vol. 103, p. 012204, Jan 2021.
  * [63] S. Pilatowsky-Cameo, D. Villaseñor, M. A. Bastarrachea-Magnani, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, “Ubiquitous quantum scarring does not prevent ergodicity,” Nature Communications, vol. 12, p. 852, Feb 2021.
  * [64] D. S. Grün, L. H. Ymai, K. Wittmann W., A. P. Tonel, A. Foerster, and J. Links, “Integrable atomtronic interferometry,” Phys. Rev. Lett., vol. 129, p. 020401, Jul 2022.
  * [65] D. S. Grün, K. Wittmann W., L. H. Ymai, J. Links, and A. Foerster, “Protocol designs for NOON states,” Commun. Phys., vol. 5, no. 1, p. 36, 2022.