---
abstract: 'We suggest to subject anharmonically trapped Bose-Einstein condensates to sinusoidal forcing with a smooth, slowly changing envelope, and to measure the coherence of the system after such pulses. In a series of measurements with successively increased maximum forcing strength one then expects an adiabatic return of the condensate to its initial state as long as the pulses remain sufficiently weak. In contrast, once the maximum driving amplitude exceeds a certain critical value there should be a drastic loss of coherence, reflecting significant heating induced by the pulse. This predicted experimental signature is traced to the loss of an effective adiabatic invariant, and to the ensuing breakdown of adiabatic motion of the system’s Floquet state when the many-body dynamics become chaotic. Our scenario is illustrated with the help of a two-site model of a forced bosonic Josephson junction, but should also hold for other, experimentally accessible configurations.'
author:
- Christoph Heinisch
- Martin Holthaus
date: 'April 27, 2016'
title: 'Entropy production within a pulsed Bose-Einstein condensate'
---

Introduction {#sec:1}
============

In 1974 an influential series of experiments on the microwave-induced multiphoton ionization of highly excited Hydrogen atoms was initiated by J. E. Bayfield and P. M. Koch [@BayfieldKoch74]. Sending a beam of fast Hydrogen atoms with principal quantum numbers ranging from $n = 63$ to $n = 69$ through an X-band microwave cavity and measuring the resulting ionization probabilities as a function of the peak amplitude of the applied oscillating electric field for a frequency of $9.9$ GHz, these authors observed the onset and saturation of ionization with increasing field strength, although, formally, the energy of about 76 microwave photons was required to reach the continuum states from the initial state $n = 66$. Remarkably, considerable ionization was observed even when the peak microwave amplitude was small by comparison with the static electric field needed to ionize the atom.

It was soon realized that the outcome of these early experiments is well described by a classical approach, based on the evolution of a representative set of classical trajectories in phase space [@LeopoldPercival78]. Since a classical anharmonic oscillator subjected to strong periodic driving naturally gives rise to chaotic dynamics, this inevitably led to the question how the chaotic behavior exhibited by the corresponding classical model is reflected in the ionization process of the real, quantum mechanical Hydrogen atom. Subsequent theoretical work on this question has led to remarkable insights. On the one hand, Casati [*et al.*]{} have predicted the existence of a certain critical microwave field strength, dubbed quantum delocalization border, above which the quantum wave packet should delocalize, and strong excitation and ionization should take place [@CasatiEtAl87]. On the other, Blümel and Smilansky have drawn attention to a qualitative change of the system’s Floquet states at the ionization border [@BlumelSmilansky87]. Significant further interest in the experiments was spurred by the hypothesis, derived from an analysis of the quantum kicked rotator [@FishmanEtAl82; @GrempelEtAl84], that the classical-quantum correspondence might be broken in the high-frequency regime: Whereas the energy of chaotic classical trajectories grows diffusively, such diffusive energy growth, and hence ionization, should be suppressed in the quantum system by means of a mechanism closely related to the Anderson localization of particles moving on a one-dimensional disordered lattice. Indeed, signatures of this high-frequency stabilization have been reported in later works [@GalvezEtAl88; @JensenEtAl91]. Moreover, a detailed investigation has been made concerning the influence of classical resonances, and their quantum mechanical counterparts, on the observed ionization signal [@KochvanLeeuwen95].

An additional twist to the interpretation of these microwave ionization experiments is provided by the fact that they are [*not*]{} performed with strictly time-periodic driving: When the fast atoms enter a microwave resonator, they experience a fringe field which, in their rest frame, corresponds to a driving amplitude which increases slowly on the time scale set by one microwave cycle [@GalvezEtAl88]. This means that the state experimented with within the cavity is actually prepared when entering it: The initial Rydberg state can either be shifted adiabatically into the connected Floquet state, or undergo multiphoton transitions at avoided crossings of its quasienergy [@BreuerHolthaus89]. In the latter case the precise form of the fringe fields, which translates into the slowly varying envelope of the pulse actually “seen” by the atoms, can strongly influence the experimental results. Therefore, the question whether or not the actual wave function of the microwave-driven Hydrogen atom is able to adiabatically follow the instantaneous Floquet states under the action of a pulse provides a key for understanding the ionization data.

In the present paper we suggest to transfer these previous single-particle experiments, which have shaped much of our current understanding of periodically driven quantum systems, to the many-body level, utilizing Bose-Einstein condensates subjected to sinusoidal driving with a smooth, slowly changing envelope. Again, the key question then is whether or not the condensate wave function responds to such pulses in an adiabatic manner. As we will argue, there may be far-reaching conceptual analogies between these older microwave experiments and the ones proposed here, with the classical-quantum correspondence being replaced by the correspondence between the mean-field description of a driven condensate and its full many-body dynamics. We proceed as follows: In Sec. \[sec:2\] we review a convenient $N$-particle model system, and discuss the results of selected numerical calculations which illustrate its interaction with forcing pulses. In Sec. \[sec:3\] we then compare the $N$-particle dynamics to the predictions of a mean-field approach, emphasizing that the latter remains trustworthy only in the regular regime, which allows for adiabatic following of the many-body wave function to the driving amplitude. In the final Sec. \[sec:4\] we formulate our conclusions.

Model calculations {#sec:2}
==================

Our numerical studies are based on the familiar model of a bosonic Josephson junction [@GatiOberthaler07] which supposes that the system consists of two sites connected by a tunneling contact, such that Bose particles sitting on the same site repel each other, while interaction between particles on different sites is neglected. Thus, the basic part of the model Hamiltonian reads [@LMG65; @ScottEilbeck86; @MilburnEtAl97; @Leggett01] $$\begin{aligned}
    H_0 & = & -\frac{\hbar\Omega}{2} 
    \left( {a^{\phantom{\dagger}}}_1{a^{\dagger}}_2 + {a^{\dagger}}_1{a^{\phantom{\dagger}}}_2 \right) 
\nonumber \\ & & 
    + \hbar\kappa \left( {a^{\dagger}}_1{a^{\dagger}}_1{a^{\phantom{\dagger}}}_1{a^{\phantom{\dagger}}}_1 
    + {a^{\dagger}}_2{a^{\dagger}}_2{a^{\phantom{\dagger}}}_2{a^{\phantom{\dagger}}}_2 \right) \; . 
\label{eq:UDJ}\end{aligned}$$ The bosonic operators ${a^{\phantom{\dagger}}}_j$ and ${a^{\dagger}}_j$, obeying the usual commutation relations $$\left[ {a^{\phantom{\dagger}}}_j, {a^{\phantom{\dagger}}}_k \right] = 0     \; , \quad 
    \left[ {a^{\dagger}}_j, {a^{\dagger}}_k \right] = 0     \; , \quad
    \left[ {a^{\phantom{\dagger}}}_j, {a^{\dagger}}_k \right] = \delta_{jk} \; ,$$ describe, respectively, the annihilation and creation of a particle at the $j$th site ($j,k = 1,2$). The tunneling matrix element is written as $\hbar\Omega/2$, so that $\Omega$ denotes the single-particle tunneling frequency, while $2\hbar\kappa$ is the repulsion energy contributed by one pair of particles occupying a common site. Although the $N$-particle ground state of this system (\[eq:UDJ\]) is exactly equal to a coherent state, that is, to an $N$-fold occupied single particle state, for vanishing interparticle interaction only, it still remains a highly coherent condensate state even for finite values of the scaled interaction strength $\alpha = N\kappa/\Omega$ [@MazzarellaEtAl11; @GertjerenkenHolthaus15b]. This is verified here by computing the one-particle reduced density matrix $$\varrho_n = \left( \begin{array}{cc}
    \langle {a^{\dagger}}_1 {a^{\phantom{\dagger}}}_1 \rangle_n & \langle {a^{\dagger}}_1 {a^{\phantom{\dagger}}}_2 \rangle_n \\
    \langle {a^{\dagger}}_2 {a^{\phantom{\dagger}}}_1 \rangle_n & \langle {a^{\dagger}}_2 {a^{\phantom{\dagger}}}_2 \rangle_n 
        \end{array} \right) \; ,$$ where the expectation values $\langle \cdots \rangle_n$ are taken with respect to the $n$th energy eigenstate, from which one obtains the invariant $$\eta_n  = 2 N^{-2} \, {\rm tr} \, \varrho_n^2 - 1 
\label{eq:ETA}$$ as a measure for the degree of coherence of that eigenstate: A pure condensate state gives ${\rm tr} \, \varrho_n^2 = N^2$, so that $\eta_n = 1$, whereas a maximally fractionalized state yields ${\rm tr} \, \varrho_n^2 = N^2/2$, amounting to $\eta_n = 0$ [@Leggett01]. In Fig. \[F\_1\] we depict this degree of coherence for the lowest five eigenstates of a junction (\[eq:UDJ\]) occupied with $N = 10\,000$ Bose particles. Evidently, the ground state $n = 0$ remains highly coherent even up to $\alpha = 2.0$, and thus serves as a good condensate state.

Many authors have previously studied the dynamics of a bosonic Josephson junction (\[eq:UDJ\]) under the action of an external time-periodic force with constant amplitude [@HolthausStenholm01; @SalmondEtAl02; @SalasnichEtAl02; @MahmudEtAl05; @WeissTeichmann08; @GertjerenkenWeiss13]. In contrast, in order to monitor the system’s adiabatic or non-adiabatic responses to forcing [*pulses*]{}, and thus to draw the parallels to the microwave ionization experiments reviewed in the Introduction, here we introduce a pulse-like, site-diagonal interaction of the form $$H_1(t) = \hbar \mu(t) \sin(\omega t)
    \left( {a^{\dagger}}_1{a^{\phantom{\dagger}}}_1 - {a^{\dagger}}_2{a^{\phantom{\dagger}}}_2 \right) \; ,$$ where $\omega$ is the carrier frequency of the pulse, and $\hbar\mu(t)$ models its envelope, assumed to be slowly varying on the scale of $T = 2\pi/\omega$. In particular, we consider Gaussian envelopes $$\mu(t) = \mu_{\rm max} \exp\left(-\frac{t^2}{2\sigma^2}\right)
\label{eq:PUL}$$ with width parameter $\sigma$; their maximum strength $\mu_{\rm max}$ thus carries the dimension of a frequency. The full model given by the total Hamiltonian $$H(t) = H_0 +  H_1(t) 
\label{eq:TOT}$$ has also been employed in further recent investigations of pulsed many-Boson dynamics [@GertjerenkenHolthaus15b; @GertjerenkenHolthaus15a; @HeinischHolthaus16]. Since the dimension of its Hilbert space figures merely as $N+1$ when the driven junction hosts $N$ Bose particles, it allows one to treat fairly large particle numbers with only moderate numerical effort.

We now stipulate that initially, at time $t/T = -\infty$, the state of the system is given by the numerically determined condensate ground state of the undriven junction (\[eq:UDJ\]) for given, fixed particle number $N$, and integrate the time-dependent Schrödinger equation with the Hamiltonian (\[eq:TOT\]) to determine the evolving $N$-particle state $| \psi(t) \rangle$. Here and in the following the reference time scale $T$ is given by the carrier cycle time $T = 2\pi/\omega$; in practice, our numerical integrations cover the interval from $t_i = -10\,\sigma$ to $t_f = +10\,\sigma$.

Figure \[F\_2\] shows the time-resolved population imbalance $$\langle J_z \rangle/N = \frac{1}{2N}
    \langle \psi(t) | {a^{\dagger}}_1{a^{\phantom{\dagger}}}_1 - {a^{\dagger}}_2{a^{\phantom{\dagger}}}_2 | \psi(t) \rangle
\label{eq:IMB}$$ of a system comprising $N = 100$ particles with scaled interaction strength $N\kappa/\Omega = 0.95$ while interacting with a pulse (\[eq:PUL\]) equipped with a carrier frequency $\omega/\Omega = 1.0$, moderate maximum strength $\mu_{\rm max}/\Omega = 0.6$, and width $\sigma/T = 5.0$. Although such a pulse actually is comparatively short, so that the variation of the envelope during one single cycle $T$ might not appear to be negligible, the system still is responding in an almost perfectly adiabatic manner: Adjusting itself to the “slowly” changing envelope, here the initial state is transformed almost perfectly into the connected instantaneous Floquet state belonging to the current driving amplitude, so that the system is able to return almost completely to its initial state with zero imbalance at the end of the pulse [@BreuerHolthaus89; @HeinischHolthaus16].

If the strength of the pulse is increased, the picture changes: In Fig. \[F\_3\] we plot the imbalance for $\mu_{\rm max}/\Omega = 0.9$; all other parameters are the same as before. Now the system’s response is no longer adiabatic: At the end of the pulse several eigenstates of the junction (\[eq:UDJ\]) are appreciably excited, leading to small, seemingly erratic fluctuations of the imbalance after the pulse is over.

In order to quantify this loss of adiabaticity with increasing pulse strength in a more systematic manner, we determine the final occupation probabilities $$p_n = \big| \langle n | \psi(t_f) \rangle \big|^2 \; ,$$ where $\{ | n \rangle \}$ denotes the eigenstates of the unperturbed system (\[eq:UDJ\]), and compute the (dimensionless) von Neumann entropy generated by the pulse according to [@PathriaBeale11] $$S = - \sum_n p_n \ln p_n \; .
\label{eq:ENT}$$ Since there are $N+1$ different states, the maximum entropy $S_{\max} = \ln(N+1) \approx \ln N$ would result if all eigenstates of the junction were populated equally after the pulse, $p_n = 1/(N+1)$ for all $n = 0, \ldots, N$. The normalized entropy $S/\ln N$ therefore varies between zero and unity, with values close to unity indicating almost equal distribution of the final state over the unperturbed energy eigenstates.

In Fig. \[F\_4\] we depict the logarithm of this normalized final entropy vs. the scaled maximum driving amplitude $\mu_{\rm max}/\Omega$, again for $N = 100$, $N\kappa/\Omega = 0.95$, and $\omega/\Omega = 1.0$. Here we also consider several different pulse widths, $\sigma/T = 5$, $10$, $20$, $30$, and $40$, thereby allowing for a varying degree of adiabaticity. Except for the shortest pulses with $\sigma/T = 5$ all curves almost coalesce, giving a quite consistent picture: Pulses with maximum scaled amplitude lower than $\mu_{\rm max}/\Omega \approx 0.65$ enable adiabatic following, and thus result in negligible entropy production. Then there is a transition regime, extending from $\mu_{\rm max}/\Omega \approx 0.65$ to about $\mu_{\rm max}/\Omega \approx 0.85$, in which a slight increase of the maximum amplitude triggers a steep rise of the entropy. For still stronger pulses one observes close-to-maximum entropy generation, corresponding to a significant heating of the initial condensate.

In Figs. \[F\_5\] and \[F\_6\] we repeat these investigations for $N = 1000$ and $N = 10\,000$ particles, respectively, keeping the scaled interaction strength $N\kappa/\Omega = 0.95$ fixed, so that an increase of the particle number $N$ is accompanied by a reduction of the bare interparticle interaction strength $\hbar\kappa$ [@GertjerenkenHolthaus15b]. Again one observes an adiabatic regime and a regime of close-to-maximum entropy production, separated by a transition regime, but there is an additional feature: With these larger particle numbers the system is able to discern the different pulse widths, meaning that the onset of entropy production shifts to higher maximum amplitudes when the pulse is made longer. On the other hand, the onset of the regime of maximum heating appears to be more or less independent of the pulses’ length.

These numerical findings bear a striking similarity to the microwave ionization experiments [@BayfieldKoch74] referred to in the Introduction: There one finds a “delocalization border” [@CasatiEtAl87] above which a microwave pulse leads to significant ionization of a highly excited Hydrogen atom; here one observes a “heating border” above which a driving pulse leads to significant entropy generation in the model defined by the Hamiltonian (\[eq:TOT\]). This model may be somewhat simple, and not capture all aspects of the exact many-body dynamics of experimentally feasible driven bosonic Josephson junctions. Nonetheless, we surmise that its main constitutive element, the existence of a heating border, is a generic feature of pulsed Bose-Einstein condensates also in typical experimental configurations which do not lend themselves to exact numerical many-body calculations.

The mean-field to $N$-particle correspondence {#sec:3}
=============================================

Within a mean field-type approach the dynamics of Bose-Einstein condensates are described by the time-dependent nonlinear Gross-Pitaevskii equation for the condensate’s macroscopic wave function [@Pitaevskii61; @Gross63; @PethickSmith08; @PitaevskiiStringari03]. In the case of the driven two-site system governed by the Hamiltonian (\[eq:TOT\]) this equation is cast into the form [@GertjerenkenHolthaus15b] $$\begin{aligned}
    {{\mathrm i}}\frac{{{\mathrm d}}}{{{\mathrm d}}\tau} c_1(\tau) & = & -\frac{1}{2} c_2(\tau)
    + 2\alpha | c_1(\tau) |^2 c_1(\tau)
\nonumber \\    & & 
    + \frac{\mu(\tau)}{\Omega}\sin\left(\frac{\omega}{\Omega}\tau\right)
    c_1(\tau) \; ,
\nonumber \\
    {{\mathrm i}}\frac{{{\mathrm d}}}{{{\mathrm d}}\tau} c_2(\tau) & = & -\frac{1}{2} c_1(\tau)
    + 2\alpha | c_2(\tau) |^2 c_2(\tau)
\nonumber \\    & &
    - \frac{\mu(\tau)}{\Omega}\sin\left(\frac{\omega}{\Omega}\tau\right)
    c_2(\tau) \; ,      
\label{eq:GPE}\end{aligned}$$ where $\tau = \Omega t$ is a dimensionless time variable, and $$\alpha = \frac{N\kappa}{\Omega}$$ denotes the scaled interparticle interaction strength already employed in the previous section. The squared amplitudes $\big| c_j(\tau) \big|^2$ are to be interpreted as the expected fraction of particles occupying the $j$th site at time $\tau$ ($j = 1,2$). A detailed derivation [@GertjerenkenHolthaus15b] of these coupled nonlinear equations (\[eq:GPE\]) which is not based on a coherent-states approach, but rather compares the evolution of an $N$-particle system to that of subsidiary systems obtained by the removal of one particle, reveals that they indeed do provide a faithful image of the exact $N$-particle dynamics generated by the Hamiltonian (\[eq:TOT\]) in the limit of large $N$, provided their solutions remain [*regular*]{}. In this case the dynamics in Fock space are stiff, in the sense that the subsidiary $(N-1)$-particle states remain closely related to the actual $N$-particle state while evolving in time. If, however, the solutions to Eqs. (\[eq:GPE\]) become [*chaotic*]{}, their direct link to the $N$-particle level is lost: In that latter case the exact $N$-particle state, when responding to the external drive, acquires a level of complexity which amounts to a destruction of the macroscopic wave function, so that its description in terms of a Gross-Pitaevskii approach becomes pointless [@GertjerenkenHolthaus15b; @GertjerenkenHolthaus15a].

We now consider the initial condition $$\left( \begin{array}{cc} 
        c_1(\tau_i) \\ c_2(\tau_i) 
    \end{array} \right) = \frac{1}{\sqrt{2}}
    \left( \begin{array}{cc}
        1 \\ 1
    \end{array} \right) \; ,$$ as corresponding to the condensate ground state $| 0 \rangle$ of the undriven junction (\[eq:UDJ\]), and utilize the Gross-Pitaevskii equation (\[eq:GPE\]) for computing the mean-field return probabilities $$P_{\rm ret}^{\rm mf} = \big|
    c_1^*(\tau_i) c_1(\tau_f) + c_2^*(\tau_i) c_2(\tau_f)
    \big|^2 
\label{eq:PMF}$$ for the same pulses as studied before. In Fig. \[F\_7\] we display such mean-field return probabilities for short and long pulses, $\sigma/T = 5$ (upper panel) and $\sigma/T = 40$ (lower panel), in comparison with the corresponding full $N$-particle quantum return probabilities $$P_{\rm ret}^{N} = \big| \langle 0 | \psi(\tau_f) \rangle \big|^2    
\label{eq:PNP}$$ obtained for $N = 10\,000$. For the shorter pulses the Gross-Pitaevskii treatment slightly overestimates the “critical” maximum amplitude $\mu_{\rm max}/\Omega$ at which the sudden drop of the $N$-particle return probability indicates the loss of adiabaticity, in line with a substantial entropy production. For the longer pulses, which should favor adiabatic motion, the onset of a chaotic mean-field return pattern agrees very well with the quantum heating border. We conclude that the comparatively simple Gross-Pitaevskii equation, although it is lacking theoretical justification in the entropy-generation regime where no macroscopic wave function exists, is quite capable of determining the heating border itself.

Discussion {#sec:4}
==========

The field of “quantum chaos”, which by now has led to a deepened understanding not only of the classical-quantum correspondence, but also of generic quantum dynamics itself, has previously concentrated on single-particle systems, chosen to be sufficiently simple to allow for full-fledged numerical analysis on both the classical and the quantum level [@Gutzwiller91; @Haake10; @Stoeckmann07]. Externally driven Bose-Einstein condensates now offer attractive prospects for extending these developments to an even more challenging many-body context. As long as such condensates respond to the drive in a “regular” manner, their macroscopic wave function remains preserved, possibly corresponding to a macroscopically occupied single-particle Floquet state. If, however, one encounters the realm of “chaotic” dynamics, the many-body state no longer remains entropyless, but necessarily has to carry a lot of information, and therefore to acquire a high level of complexity which is not accessible to a mean-field ansatz.

The mean-field dynamics associated with a periodically driven bosonic two-mode Josephson junction (\[eq:UDJ\]) are exactly equivalent to those of a driven nonlinear classical pendulum [@GertjerenkenHolthaus15b]. Therefore, in this exceptional case one may approach chaotic many-body dynamics by mapping out the corresponding classical phase space, as done in a recent experiment reported by Tomkovič [*et al.*]{} [@TomkovicEtAl15]. However, in more general cases a reduction to two-mode dynamics will not be feasible, and more versatile tools of investigation will be required. Then the concept of adiabatic following of driven many-body states to the respective instantaneous Floquet states, as explored in the context of “avoided-level-crossing spectroscopy” for driven Bose-Hubbard-like systems [@EckardtHolthaus08], suggests itself: Adiabatic following of a Bose-Einstein condensate to a smooth forcing pulse amounts to the preservation of an entropyless condensate state which can exist in the regular regime only, whereas the entropy production associated with the loss of adiabaticity signals the onset of chaos.

Thus, it should be a worthwhile enterprise to set up an experiment for investigating the transition from regular, mean field-like condensate dynamics to chaotic many-body dynamics in systematic detail along the lines suggested in the present proposal: Start from a Bose-Einstein condensate in a trapping configuration which allows for the application of well-defined forcing pulses. Then subject the condensate to a pulse with a slowly varying, smooth envelope, and measure the coherence of the final state after the pulse by recording, [*e.g.*]{}, the condensate fraction through time-of-flight absorption imaging. As long as the maximum pulse strength does not reach the chaotic regime, the system’s response should be adiabatic. Thus, the condensate will return almost fully to its initial state after the pulse, with negligible production of entropy. If, however, the driving amplitude crosses the border to chaos, adiabatic following is disabled, and the system undergoes manifold excitations during the pulse, corresponding to a high entropy production. Such dynamical heating of the condensate leads to its destruction, detectable after the pulse through the loss of coherence. Therefore, in configurations possessing a sharp chaos border, which may be achievable with anharmonic trapping potentials, one should observe a sharp drop of the final coherence when the maximum driving amplitude is gradually increased.

From a theoretician’s perspective, the difference between “regular” and “chaotic” dynamics of periodically driven quantum systems is reflected in the properties of their quasienergy spectra when these are considered for all instantaneous driving amplitudes encountered during the pulses, with smooth coarse-grained quasienergy eigenvalues in the regular regime enabling effectively adiabatic following of the corresponding Floquet states, whereas multiple avoided quasienergy crossings in the chaotic regime give rise to multiphoton-like Landau-Zener transitions effectuating the observed entropy growth [@Holthaus16]. While the numerical computation of such quasienergy spectra, and the verification of adiabatic following of many-body Floquet states, is still feasible for simplified models such as the driven Josephson junction considered here [@HeinischHolthaus16], this goal will remain unachievable for realistic mesoscopic condensates consisting of $10^6$ particles, say. Hence, the experiments suggested in this work might break new ground in an area where numerical guidance is feasible on the mean-field level only: The question whether the educated guesses derived from the present model study actually hold water for very large $N$ has to be decided in the laboratory.

Such condensate experiments would constitute a natural extension of previous microwave ionization experiments with highly excited Hydrogen atoms [@BayfieldKoch74], in which a classical chaos border manifests itself through the onset of strong ionization. While that ionization border would find its match in the heating border discussed here, the question whether driven condensates may also give rise to an analog of the Anderson-like suppression of classical phase-space diffusion reported for the quantum kicked rotator [@FishmanEtAl82; @GrempelEtAl84] appears to be open.

We are grateful for CPU time granted to us on the HPC cluster HERO, located at the University of Oldenburg and funded by the DFG through its Major Research Instrumentation Programme (INST 184/108-1 FUGG), and by the Ministry of Science and Culture (MWK) of the Lower Saxony State.

[99]{}

J. E. Bayfield and P. M. Koch, Phys. Rev. Lett. [**33**]{}, 258 (1974).

J. G. Leopold and I. C. Percival, Phys. Rev. Lett. [**41**]{}, 944 (1978).

G. Casati, B. V. Chirikov, D. L. Shepelyansky, and I. Guarneri, Phys. Rep. [**154**]{}, 77 (1987).

R. Blümel and U. Smilansky, Z. Phys. D [**6**]{}, 83 (1987).

S. Fishman, D. R. Grempel, and R. E. Prange, Phys. Rev. Lett. [**49**]{}, 509 (1982)

D. R. Grempel, R. E. Prange, and S. Fishman, Phys. Rev. A [**29**]{}, 1639 (1984).

E. J. Galvez, B. E. Sauer, L. Moorman, P. M. Koch, and D. Richards, Phys. Rev. Lett. [**61**]{}, 2011 (1988).

R. V. Jensen, S. M. Susskind, and M. M. Sanders, Phys. Rep. [**201**]{}, 1 (1991).

P.M. Koch and K.A.H. van Leeuwen, Phys. Rep. [**255**]{}, 289 (1995).

H. P. Breuer and M. Holthaus, Z. Phys. D [**11**]{}, 1 (1989).

R. Gati and M. K. Oberthaler, J. Phys. B: At. Mol. Opt. Phys. [**40**]{}, R61 (2007).

H. J. Lipkin, N. Meshkov, and A. J. Glick, Nucl. Phys. [**62**]{}, 188 (1965).

A. C. Scott and J. C. Eilbeck, Phys. Lett. A [**119**]{}, 60 (1986).

G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, Phys. Rev. A [**55**]{}, 4318 (1997).

A. J. Leggett, Rev. Mod. Phys. [**73**]{}, 307 (2001).

G. Mazzarella, L. Salasnich, A. Parola, and F. Toigo, Phys. Rev. A [**83**]{}, 053607 (2011).

B. Gertjerenken and M. Holthaus, Annals of Physics [**362**]{}, 482 (2015).

M. Holthaus and S. Stenholm, Eur. Phys. J. B [**20**]{}, 451 (2001).

G. L. Salmond, C. A. Holmes, and G. J. Milburn, Phys. Rev. A [**65**]{}, 033623 (2002).

L. Salasnich, A. Parola, and L. Reatto, J. Phys. B: At. Mol. Opt. Phys. [**35**]{}, 3205 (2002).

K. W. Mahmud, H. Perry, and W. P. Reinhardt, Phys. Rev. A [**71**]{}, 023615 (2005).

C. Weiss and N. Teichmann, Phys. Rev. Lett. [**100**]{}, 140408 (2008).

B. Gertjerenken and C. Weiss, Phys. Rev. A [**88**]{}, 033608 (2013).

B. Gertjerenken and M. Holthaus, EPL [**111**]{}, 30006 (2015).

C. Heinisch and M. Holthaus, J. Mod. Opt. (2016); DOI: 10.1080/09500340.2016.1167263.

See, e.g., R. K. Pathria and Paul D. Beale, [*Statistical Mechanics*]{}. Third Edition (Elsevier / Butterworth-Heinemann, Oxford, 2011).

L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. [**40**]{}, 646 \[Sov. Phys. JETP [**13**]{}, 451\] (1961).

E. P. Gross, J. Math. Phys. [**4**]{}, 195 (1963).

C. J. Pethick and H. Smith, [*Bose-Einstein Condensation in Dilute Gases*]{}. Second Edition (Cambridge University Press, Cambridge, 2008).

L. Pitaevskii and S. Stringari, [*Bose-Einstein Condensation*]{} (Clarendon Press, Oxford, 2003).

M. C. Gutzwiller, [*Chaos in Classical and Quantum Mechanics*]{} (Springer-Verlag, New York Berlin Heidelberg, 1991).

F. Haake, [*Quantum Signatures of Chaos*]{}. Third Edition (Springer-Verlag, Berlin Heidelberg, 2010).

H.-J. Stöckmann, [*Quantum Chaos – An Introduction*]{} (Cambridge University Press, Cambridge, 1999).

J. Tomkovič, W. Muessel, H. Strobel, S. Löck, P. Schlagheck, R. Ketzmerick, and M. K. Oberthaler, arXiv:1509.01809.

A. Eckardt and M. Holthaus, Phys. Rev. Lett. [**101**]{}, 245302 (2008).

For a tutorial review, see: M. Holthaus, J. Phys. B: At. Mol. Opt. Phys. [**49**]{}, 013001 (2016).
