Classical and quantum chaos in a three-mode bosonic system

TL;DR

This paper explores the connection between classical chaos and quantum signatures in a three-mode bosonic system, demonstrating how classical dynamics influence quantum level statistics and out-of-time-ordered correlators, with potential experimental realizations.

Contribution

It introduces a tunable three-mode bosonic model linking classical chaos with quantum signatures, including level statistics and OTOCs, and proposes an experimental setup with Bose-Einstein condensates.

Findings

Classical chaos affects quantum level statistics, changing from Poissonian to Wigner-Dyson.

OTOCs differ significantly between regular and chaotic classical regions.

Semi-classical simulations accurately capture short-time quantum dynamics.

Abstract

We study the dynamics of a three-mode bosonic system with mode-changing interactions. For large mode occupations the short-time dynamics is well described by classical mean-field equations allowing us to study chaotic dynamics in the classical system and its signatures in the corresponding quantum dynamics. By introducing a symmetry-breaking term we tune the classical dynamics from integrable to strongly chaotic which we demonstrate by calculating Poincar\'e sections and Lyapunov exponents. The corresponding quantum system features level statistics that change from Poissonian in the integrable to Wigner-Dyson in the chaotic case. We investigate the behavior of out-of-time-ordered correlators (OTOCs), specifically the squared commutator, for initial states located in regular and chaotic regions of the classical mixed phase space and find marked differences between the two cases. The short-time behavior is well captured by semi-classical truncated Wigner simulations directly relating these features to properties of the underlying classical mean field dynamics. We discuss a possible experimental realization of this model system in a Bose-Einstein condensate of rubidium atoms which allows reversing the sign of the Hamiltonian required for measuring OTOCs experimentally.

Equations

$$C(t)\to {\left|\left\{x(t),p(0)\right\}\right|}^2={\left|\frac{\partial x(t)}{\partial x(0)}\right|}^2\sim e^{2\lambda t}.$$

$$\begin{array}{rl}\hat{H}& =g\{({\hat{a}}_0^{†}{\hat{a}}_0^{†}{\hat{a}}_1{\hat{a}}_{-1}+{\hat{a}}_1^{†}{\hat{a}}_{-1}^{†}{\hat{a}}_0{\hat{a}}_0)+{\hat{N}}_0({\hat{N}}_1+{\hat{N}}_{-1})\\ & +\frac{1}{2}({\hat{N}}_1-{\hat{N}}_{-1}{)}^2\}+q({\hat{N}}_1+{\hat{N}}_{-1})\\ & +\frac{r}{\sqrt{2}}\left\{({\hat{a}}_1^{†}+{\hat{a}}_{-1}^{†}){\hat{a}}_0+{\hat{a}}_0^{†}({\hat{a}}_1+{\hat{a}}_{-1})\right\}\end{array}$$

$$\begin{array}{rl}{H}_{\text{mf}}& =gN\{{\zeta }_0^{\ast 2}{\zeta }_1{\zeta }_{-1}+{\zeta }_1^{\ast }{\zeta }_{-1}^{\ast }{\zeta }_0^2+|{\zeta }_0{|}^2\left(|{\zeta }_1{|}^2+|{\zeta }_{-1}{|}^2\right)\\ & +m^2/2\}+q\left(|{\zeta }_1{|}^2+|{\zeta }_{-1}{|}^2\right)\\ & +r/\sqrt{2}\left\{({\zeta }_1^{\ast }+{\zeta }_{-1}^{\ast }){\zeta }_0+{\zeta }_0^{\ast }({\zeta }_1+{\zeta }_{-1})\right\}.\end{array}$$

$$\begin{array}{rl}i{\dot{\zeta }}_1& =q{\zeta }_1+gN\left\{({\rho }_1+{\rho }_0-{\rho }_{-1}){\zeta }_1+{\zeta }_{-1}^{\star }{\zeta }_0^2\right\}+\frac{r}{\sqrt{2}}{\zeta }_0\\ i{\dot{\zeta }}_0& =gN\left\{({\rho }_1+{\rho }_{-1}){\zeta }_0+2{\zeta }_0^{\star }{\zeta }_1{\zeta }_{-1}\right\}+\frac{r}{\sqrt{2}}({\zeta }_1+{\zeta }_{-1})\\ i{\dot{\zeta }}_{-1}& =q{\zeta }_{-1}+gN\left\{(-{\rho }_1+{\rho }_0+{\rho }_{-1}){\zeta }_{-1}+{\zeta }_1^{\star }{\zeta }_0^2\right\}\\ & +\frac{r}{\sqrt{2}}{\zeta }_0\end{array}$$

$$\begin{array}{rl}{\zeta }_1& =\sqrt{\frac{1-{\rho }_0+m}{2}}{e}^{i\frac{{\Theta }_s+{\Theta }_m}{2}}\\ {\zeta }_0& =\sqrt{{\rho }_0}\\ {\zeta }_{-1}& =\sqrt{\frac{1-{\rho }_0-m}{2}}{e}^{i\frac{{\Theta }_s-{\Theta }_m}{2}},\end{array}$$

$$\begin{array}{rl}{H}_{\text{mf}}& =gN{\rho }_0\left\{\left(1-{\rho }_0\right)+\sqrt{{\left(1-{\rho }_0\right)}^2-m^2}\cos {\Theta }_s\right\}\\ & +\frac{gN}{2}{m}^2+q\left(1-{\rho }_0\right)\\ & +r\sqrt{{\rho }_0}\{\sqrt{1-{\rho }_0+m}\cos \left(\frac{{\Theta }_s+{\Theta }_m}{2}\right)\\ & +\sqrt{1-{\rho }_0-m}\cos \left(\frac{{\Theta }_s-{\Theta }_m}{2}\right)\}.\end{array}$$

$$\frac{\left|\Delta \vec{x}(t)\right|}{\left|\Delta {\vec{x}}_0\right|}\sim e^{\lambda t}.$$

$$\lambda =\underset{t\to \infty }{\lim }\underset{|\Delta {\vec{x}}_0|\to 0}{\lim }\frac{1}{t}\log \frac{\left|\Delta \vec{x}(t)\right|}{\left|\Delta {\vec{x}}_0\right|}.$$

$$\lambda =\underset{t\to \infty }{\lim }\frac{1}{t}\log \frac{\left|\vec{\xi }(t)\right|}{\left|\vec{\xi }(0)\right|}.$$

$$\left(\begin{array}{c}\dot{\vec{\zeta }}\\ \dot{\Phi }\end{array}\right)=\left(\begin{array}{c}\vec{F}(\vec{\zeta })\\ D_{\vec{\zeta }}\vec{F}(\vec{\zeta })\Phi (\vec{\zeta })\end{array}\right).$$

$$P_b(s)=\alpha (b+1)s^b\exp \left[-\alpha s^{b+1}\right],\alpha ={\left[\Gamma \left(\frac{b+2}{b+1}\right)\right]}^{b+1}.$$

$$Q(\alpha )={\left|⟨\alpha |\psi ⟩\right|}^2$$

$$⟨{\left|\left[{\hat{\rho }}_0(t),{\hat{\rho }}_0(0)\right]\right|}^2⟩={\hslash }_{\text{eff}}^2{⟨{\left|2{\zeta }_0^R(0)\frac{\partial {\rho }_0(t)}{\partial {\zeta }_0^I}\right|}^2⟩}_W$$

$$\begin{array}{rl}⟨\hat{W}(t{)}^{†}\hat{V}\hat{W}(t)⟩& =⟨\hat{V}⟩+i\varphi ⟨[\hat{A}(t),\hat{V}]⟩\\ & +{\varphi }^2\left(⟨\hat{A}(t)\hat{V}\hat{A}(t)⟩-\frac{1}{2}⟨\{\hat{A}(t{)}^2,\hat{V}\}⟩\right)\\ & +\mathcal{O}({\varphi }^3)\end{array}$$

$$\begin{array}{rl}C(t)& =⟨[\hat{A}(t),\hat{V}{]}^{†}[\hat{A}(t),\hat{V}]⟩\\ & =-2\Lambda ⟨\hat{A}(t)\hat{V}\hat{A}(t)⟩+⟨\hat{A}(t){\hat{V}}^2\hat{A}(t)⟩+{\Lambda }^2⟨{\hat{A}}^2(t)⟩\\ & =⟨\hat{A}(t)(\hat{V}-\Lambda {)}^2\hat{A}(t)⟩.\end{array}$$

$$\begin{array}{rl}⟨\hat{W}(t{)}^{†}{\hat{V}}^n\hat{W}(t)⟩& ={\Lambda }^n+{\varphi }^2\left[⟨\hat{A}(t){\hat{V}}^n\hat{A}(t)⟩-{\Lambda }^n⟨{\hat{A}}^2(t)⟩\right]\\ & \equiv {\Lambda }^n+{\varphi }^2{\Gamma }_{{\hat{V}}^n}^{(2)}\end{array}$$

$$C(t)=-2\Lambda {\Gamma }_{\hat{V}}^{(2)}+{\Gamma }_{{\hat{V}}^2}^{(2)}$$

$$\begin{array}{rl}{\left[\hat{A},\hat{B}\right]}_W& =i{\hslash }_{\text{eff}}{\left\{A_W,B_W\right\}}_{MB}\\ & =i{\hslash }_{\text{eff}}\left\{A_W,B_W\right\}+\mathcal{O}\left({\hslash }_{\text{eff}}^2\right)\end{array}$$

$$\left\{A,B\right\}\equiv \sum _{i=1,0,-1}\frac{\partial A}{\partial {\zeta }_i^R}\frac{\partial B}{\partial {\zeta }_i^I}-\frac{\partial A}{\partial {\zeta }_i^I}\frac{\partial B}{\partial {\zeta }_i^R}$$

$$\begin{array}{rl}& ⟨{\left|\left[\widehat{{\rho }_0}(t),\widehat{{\rho }_0}(0)\right]\right|}^2⟩={\hslash }_{\text{eff}}^2{⟨{\left|\left\{{\rho }_0(t),{\rho }_0(0)\right\}\right|}^2⟩}_W=\\ & {\hslash }_{\text{eff}}^2{⟨{\left|\frac{\partial {\rho }_0(t)}{\partial {\zeta }_0^I}\frac{\partial {\rho }_0(0)}{\partial {\zeta }_0^R}\right|}^2⟩}_W={\hslash }_{\text{eff}}^2{⟨{\left|2{\zeta }_0^R(0)\frac{\partial {\rho }_0(t)}{\partial {\zeta }_0^I}\right|}^2⟩}_W\end{array}$$

$$\frac{\partial {\rho }_0(t)}{\partial {\zeta }_0^I}=2\left\{{\zeta }_0^R(t)\frac{\partial {\zeta }_0^R(t)}{\partial {\zeta }_0^I}+{\zeta }_0^I(t)\frac{{\zeta }_0^I(t)}{\partial {\zeta }_0^I}\right\}$$

Full text

I Introduction

Statistical mechanics is based on the assumption that in thermal equilibrium all energetically accessible microstates of a system are populated with equal probability. This assumption allows to derive a comprehensive description of equilibrium phenomena. In contrast, the question how equilibrium is reached if the system is prepared in a state far from equilibrium turns out to be rather difficult to answer in general. For classical systems the answer is closely connected to the notion of classical chaos. If trajectories in phase space are not closed and a given initial state eventually explores the whole accessible phase space volume, then the system thermalizes in the sense that the long-time average of an observable is equal to the microcanonical average Penrose (1979). This notion of chaos also entails the exponential sensitivity to the initial condition, i.e. that initially close trajectories will deviate exponentially from each other at later times, which can be quantified by the Lyapunov exponent Eckmann and Ruelle (2004).

For quantum systems the relation between the microscopic dynamics and thermodynamic ensembles is more subtle since the notion of trajectories in phase space is in general absent due to Heisenberg’s uncertainty principle. Moreover, the dynamics governed by Schrödinger’s equation actually conserves the overlap between two initial states. Thus, a crucial step towards understanding how and in what sense closed quantum systems thermalize is to extend the notion of chaos and ergodicity to the quantum domain in a consistent way Haake (2010); D’Alessio et al. (2016). For systems which feature a well-defined classical limit it has been shown that the statistics of the energy levels of the quantum Hamiltonian are a faithful indicator of chaos in the corresponding classical model. Wigner and Dyson famously noticed that for a generic quantum many-body system the distribution of energy levels in a small energy window essentially looks like that of a random Hamiltonian matrix Wigner (1993); Dyson (1962). This insight led to the conjecture that the level statistics of quantum systems that have a classically chaotic counterpart can be described by random matrix theory Bohigas et al. (1984) which has been confirmed numerically Brody et al. (1981); Rudnick (2008); Wintgen and Friedrich (1987); Rigol and Santos (2010); Santos et al. (2012).

Despite these successes, the statistics of the eigenstates of a quantum Hamiltonian only provides a static picture, which limits its usefulness for answering questions about relaxation dynamics in quantum systems. Instead, the spreading of correlations and entanglement Calabrese et al. (2009); Calabrese and Cardy (2005) has shifted into the focus as a tool for understanding relaxation in closed quantum systems. Recently, a specific type of correlation functions, so-called out-of-time-order correlators (OTOCs), have been proposed as suitable indicators of chaos in quantum systems Cotler et al. (2018); Rammensee et al. (2018); Larkin and Ovchinnikov (1969); Schuckert and Knap (2019). In particular the squared commutator $C(t)=\langle[\hat{W}(t),\hat{V}(0)]^{\dagger}[\hat{W}(t),\hat{V}(0)]\rangle$ 111Note that sometimes only the out-of-time order part $\hat{W}(t)\hat{V}(0)\hat{W}(t)\hat{V}(0)$ is defined as OTOC, however, we will use the term OTOC also when referring to the squared commutator. is of interest as its growth can be related to the Lyapunov exponent by heuristic semi-classical arguments. Here, $\hat{V}(0)$ and $\hat{W}(t)=\hat{U}(t)^{\dagger}\hat{W}\hat{U}(t)$ are Heisenberg operators 222The name OTOC comes from the fact that upon expanding the squared commutator one obtains correlation functions in which the Heisenberg operators do not act in normal time order, i.e. cannot be put on a single Keldysh contour.. For example, if $\hat{V}$ and $\hat{W}$ correspond to position and momentum operators $\hat{x}$ and $\hat{p}$, in the limit of small $\hbar$ the short-time behavior of $C(t)$ is obtained by replacing operators by classical fields and the commutator by the Poisson bracket (see Refs. Cotler et al. (2018) for a more rigorous argument). This results in the quantum to classical correspondence Cotler et al. (2018); Rammensee et al. (2018); Rozenbaum et al. (2017); Chávez-Carlos et al. (2019)

[TABLE]

where $\lambda$ is the classical Lyapunov exponent. This simple argument suggests that the squared commutator grows exponentially at short times up to the Ehrenfest time Rammensee et al. (2018) if the corresponding classical model is chaotic. For this reason the growth rate of the squared commutator has been termed the quantum Lyapunov exponent. We emphasize that the mentioned semi-classical arguments only hold in systems with a well-defined classical limit of high occupation numbers. Also, as we discuss in more detail below, in higher dimensional systems the derivative appearing in Eq. (1) will, for a given OTOC, involve a specific pair of phase space coordinates [$\partial x_{i}(t)/\partial x_{j}(0)$] and thus will not necessarily grow exponentially with the largest Lyapunov exponent of the classical dynamics, which characterizes of the phase space direction of fastest growth. OTOCs have also been of interest in the context of spin systems with local interactions and random unitary circuits where they quantify operator spreading, or scrambling of quantum information. The present study focuses on the case of a bosonic few mode system where a classical mean-field description becomes exact is the limit of large particle numbers $N$.

Intense theoretical efforts have recently been undertaken to reach a more thorough understanding of the squared commutator as an indicator for quantum chaos (see Swingle (2018) for a recent review). Measuring OTOCs experimentally has been challenging due to the requirement of implementing many-body echo protocols Gärttner et al. (2017); Li et al. (2017); Wei et al. (2018); Niknam et al. (2020). In particular, the predicted exponential growth of OTOCs in systems showing classical chaos in their mean-field dynamics has not been observed in experiments. This motivates us to study OTOCs in a model of three bosonic modes which can be realized experimentally by a Bose-Einstein condensate (BEC) of rubidium atoms. We demonstrate that this model shows rich dynamics being tunable between regimes of regular and chaotic dynamics and offers a way to implement the many-body echo experimentally, thus providing an ideal platform for studying OTOCs and quantum chaos.

Specifically, we consider a BEC in a tightly confining trapping potential which constrains the dynamics to the lowest trap state. The relevant degrees of freedom are the three Zeeman levels of the $F=1$ hyperfine manifold of 87Rb. The dynamics of this three-mode system is governed by mode-changing collisions between the atoms and the quadratic Zeeman shift. The classical mean field dynamics of this model is integrable and has been studied extensively including several experimental realizations Huang et al. (2015); Plimak et al. (2006); Law et al. (1998); Kronjäger et al. (2005); Schmaljohann et al. (2004); Hamley et al. (2012); Lange et al. (2018); Kunkel et al. (2018); Linnemann et al. (2016); Gerving et al. (2012); Davis et al. (2019). By introducing a tunable coherent coupling term between neighboring Zeeman states Kunkel et al. (2019) integrability is broken 333Recall that we refer to the integrability of classical mean field equations. Integrability is broken by a symmetry-breaking field. This is different from the often considered case of an integrable system of non-interacting particles where integrability is broken by adding interactions between the particles..

The objective of this work is thus to provide a detailed understanding of the classical mean field dynamics of this model and to explore signatures of classical chaos in the quantum dynamics at finite particle number $N$. For this we first characterize the classical dynamics by calculating Poincaré sections and Lyapunov exponents for varying strength of the integrability breaking term. We find that the classical phase space portraits and Lyapunov exponents show a transition from regular to mixed phase space and strongly chaotic dynamics. Turning to quantifiers of quantum chaotic behavior we find that the level statistics changes from Poissonian to Wigner-Dyson-like (more precisely to a Brody parameter of $\approx 0.6$ Brody et al. (1981)). Finally, the squared commutator is found to generically show a faster growth for initial states located in classically chaotic regions than in regular ones. The short-time behavior before the Ehrenfest time, or scrambling time, $t_{s}\propto\log(N)$, is well described by semi-classical methods showing that the OTOC indeed reveals signatures of classical chaos in the quantum dynamics at finite $N$. We find that OTOCs with respect to an initial coherent state located in a classically chaotic region of phase space often show exponential growth at short times, however, the precise behavior is expected to depend strongly on the initial state and will, in general, not yield the largest Lyapunov exponent of the limiting classical dynamics.

Related work on OTOCs and chaos in few mode bosonic systems includes a study reporting on the growth of the squared commutator in the integrable limit of a model closely related to the one considered here Hummel et al. (2019). The squared commutator is found to show exponential growth if the system is initially prepared at a dynamically unstable point of the classical phase space. Classically chaotic dynamics has been found to emerge when extending this model to five modes Kronjäger et al. (2008), when considering the motional modes of a single component BEC in an anharmonic trapping potential Garcia-March et al. (2018), and in a three-site Bose-Hubbard model (Bose trimer) Franzosi and Penna (2003); Nemoto et al. (2000); Buonsante et al. (2009). In the latter cases the interactions are quite different from the mode-changing collisions considered here, leading to different types of phase space structures. OTOCs have not yet been discussed in the case of few mode systems with classically chaotic dynamics.

The remainder of this work is structured as follows: After introducing the model Hamiltonian in Sec. II, we first analyse the classical phase space to determine the parameter regimes in which chaotic behaviour is encountered (Sec. III). Then, in Sec. IV we conduct the corresponding analysis for the quantum case, looking for quantum signatures of classical chaos. Section V details the prospects of measuring OTOCs experimentally in a spinor BEC. We draw our conclusions and discuss possible future research directions in Sec. VI.

II Model

We consider the following three-mode Hamiltonian in second quantized form Gerving et al. (2012)

[TABLE]

where $\hat{a}_{i}$ ($\hat{a}_{i}^{\dagger}$), with $i=0,\pm 1$, denote annihilation (creation) operators of three bosonic modes giving rise to the corresponding number operators $\hat{N}_{i}\equiv\hat{a}_{i}^{\dagger}\hat{a}_{i}$. In the experimental realization of this model in a spinor BEC Stamper-Kurn and Ueda (2013); Gerving et al. (2012); Hamley et al. (2012); Lange et al. (2018); Kunkel et al. (2018); Linnemann et al. (2016) mentioned in Sec. I the parameter $g$ represents the strength of s-wave interactions between the atoms, $q$ is the quadratic Zeeman shift, and $r$ corresponds to a coherent coupling between the modes induced by an rf-modulation of the magnetic field. All these parameters are tunable. For more details on the proposed experimental realization see Sec. V.

The Hamiltonian (2), in addition to the total energy, conserves the total number of particles $\hat{N}\equiv\hat{N}_{-1}+\hat{N}_{0}+\hat{N}_{1}$ and for $r=0$ it also conserves the imbalance, or magnetization, $\hat{S}_{z}\equiv\hat{N}_{1}-\hat{N}_{-1}$. Thus, as we will see shortly, at $r=0$ the model is integrable in the classical mean field limit. By tuning $r\neq 0$ magnetization is no longer conserved and integrability is broken. In the following we fix $gN$ as the reference energy scale. Since we also set $\hbar=1$, times are given in units of $\hbar/gN$. The observed qualitative features are independent of the exact choice of $q/gN$ as long as it is at most of order $1$. We thus fix $gN=q=1$ in the following discussion. The parameter $r$ will be varied from the classically integrable case ($r=0$) to the strongly chaotic case ($r\approx 1$).

Bosonic systems with a finite number of modes can be described by classical mean-field equations in the limit of large mode occupations and at times shorter than the Ehrenfest time Sciolla and Biroli (2011), which scales as $\sqrt{N}$ for integrable and $\log(N)$ for classically chaotic systems. For finite $N$ the dynamics can be described as that of classical wave packets of width $\propto 1/\sqrt{N}$ and thus the classical limit is controlled by $\hbar_{\text{eff}}\sim 1/N$ Pappalardi et al. (2018). Technically, the classical mean field Hamiltonian is obtained by replacing the creation and annihilation operators $\hat{a}_{i}$ and $\hat{a}^{\dagger}_{i}$ by $\zeta_{i}\sqrt{N}$ and $\zeta_{i}^{\star}\sqrt{N}$, respectively, where $\zeta_{i}$ are complex-valued, normalized classical fields. This approximation becomes exact in the large-$N$ limit where terms of order $1/N$ are negigible Gardiner and Zoller (2004). The resulting mean-field Hamiltonian reads (keeping $gN$ fixed)

[TABLE]

From Eq. (3) one straightforwardly derives the equations of motion (setting $\hbar=1$)

[TABLE]

where $\rho_{i}=|\zeta_{i}|^{2}$. We integrate these equations numerically in order to find the classical time evolution.

It is convenient to reparametrize the fields $\zeta_{i}$ by explicitly exploiting conservation laws Gerving et al. (2012): We express the fields by real amplitudes and phases $\zeta_{j}=\sqrt{\rho_{j}}\ e^{i\varTheta_{j}}$. As mentioned above, the Hamiltonian conserves the total particle number $N$ which implies that the norm of $\sum_{i}\rho_{i}$ is constant ($=1$). Fixing also the global phase reduces the number of relevant classical variables (phase space coordinates) from six to four Hamley et al. (2012)

[TABLE]

with the normalized magnetization $m=\rho_{1}-\rho_{-1}=(N_{1}-N_{-1})/N$, the Larmor precession phase $\varTheta_{m}=\varTheta_{1}-\varTheta_{-1}$ and the spinor phase $\varTheta_{s}=\varTheta_{1}+\varTheta_{-1}-2\varTheta_{0}$. As we will use this parametrization in the following description of the classical phase space, it is useful to rewrite the mean field Hamiltonian as

[TABLE]

As the canonical transformation (5) leaves us with two pairs of canonical coordinates $\left\{\rho_{0},\varTheta_{S}\right\}$ and $\left\{m,\varTheta_{m}\right\}$, we immediately see from Eq. (6) that the magnetization $m$ is conserved for $r=0$.

III Classical phase space

In this section we characterize the classical dynamics and study its dependence on the strength of the integrability breaking term. We identify regimes of classically chaotic dynamics adopting the following notion of classical chaos. For a system to be chaotic we require that it is both sensitive to initial conditions, meaning that two initially close-by trajectories diverge exponentially, and also that any two trajectories with equal energy and possibly other conserved quantities come arbitrarily close to each other at later times. While the first requirement can be checked numerically by calculating Lyapunov exponents, the second is implied by ergodicity (almost every trajectory passes arbitrarily close to almost every point in phase space) Meredith et al. (1988) which can be numerically verified by studying Poincaré sections.

III.1 Poincaré sections

Instead of looking at the full four dimensional trajectories in classical phase space, we only record their crossings with a certain (three dimensional) surface of section, here defined by $\varTheta_{m}=0$. For visualization, we further use that in a time-independent classical Hamiltonian system energy is conserved, such that we can project our Poincaré section onto a two dimensional plane (in this case the $\left\{\rho_{0},\varTheta_{s}\right\}$-plane, see Fig. 1).

The Poincaré sections contain all relevant information to determine the range of parameters and initial states in which the system exhibits ergodic behaviour Berry (1978). When the intersection points are strictly confined to a one-dimensional subspace of the section, dynamics is regular. When the points of intersection are instead spread over a larger area, only limited by conservation of energy and norm, the system is likely to be ergodic in the sense described above. In a chaotic system the described phase space properties can be found for arbitrary surfaces of section but it turns out that the choice of the hyperplane influences how well the transition from regular to chaotic behaviour can be visualized. We found the surface defined by $\varTheta_{m}=0$ to be a suitable choice in the case at hand.

Figure 1(a) shows a sample of trajectories for the regular case $r=0$. In this case $m$ is conserved and thus $\theta_{m}$ is cyclic. To be able to show multiple Poincaré sections in one plot, each initial state (colour encoded) corresponds to a different energy. The Poincaré sections in this regular case are confined to one dimension (lines). Figures 1(b) and (c) show the Poincaré sections for non-integrable cases. We used a raster of initial states in the $\left\{\rho_{0},\varTheta_{s}\right\}$-plane with $m$ and $\varTheta_{m}$ chosen such that all trajectories have the same fixed energy of $E=1.005$. This energy shell was chosen in order to ensure that the energetically allowed region in the $\left\{\rho_{0},\varTheta_{s}\right\}$-plane stays sufficiently large for all considered values of $r$. It matches the energy of the separatrix between vibrational and rotational motion [solid black line in Fig. 1(a)] in the integrable case ($r=0$).

We observe the transition from a regular [Fig. 1(a)] to a strongly chaotic phase space [Fig. 1(c)], transitioning a mixed phase space [Fig. 1(b)] where depending on the initial state, a trajectory is confined to a one dimensional subspace or samples a larger part of the available phase space. The regular structure first starts to dissolve around the separatrix. The regions around the stable fixed points of the integrable dynamics retain closed orbits up to relatively large values of $r$. In the case of $r=0.5$, where we do not observe any regular islands any more in the shown energy shell, we still found mixed phase space structure at other energies, which is why we do not call this case fully chaotic but rather strongly chaotic.

We note here that we computed Poincaré sections for a large set of values of $q$ and $r$ (and accordingly adjusted energy $E$) in addition to the ones shown in Fig. 1 and found that the qualitative features are the same in all cases: If one of the parameters dominates (e.g. $q\gg 1$ or $r\gg 1$), or if $r\approx 0$, the system looks regular. If all parameters are of similar magnitude, ergodic parts of the phase space emerge and eventually dominate. This transition is not as well visible for all choices of parameters and energy since regions where regular and chaotic patches coexist may only be found in certain parts of the phase space.

III.2 Lyapunov exponents

As an ergodic phase space is only a necessary condition for classical chaos, we also need to check for the linear instability of the trajectories, measured by the Lyapunov exponent. The sensitivity of a dynamical system to initial conditions can be quantified by the (largest) Lyapunov exponent $\lambda$ Cotler et al. (2018). Given a classical trajectory in phase space $\vec{x}(t)$ and some initial conditions $\vec{x}_{0}$, the sensitivity to initial conditions can be expressed as the exponential divergence of initially infinitesimally separated trajectories:

[TABLE]

Thus, the largest Lyapunov exponent $\lambda$ is defined as

[TABLE]

For chaotic systems the largest Lyapunov exponent $\lambda$ is greater than zero, i.e. we are interested in the long-time limit as this will converge to the largest Lyapunov exponent for almost all orientations of the initial separation Eckmann and Ruelle (2004).

For the numerical calculations we use an equivalent definition of the Lyapunov exponent explicitly in terms of two (infinitesimally) close trajectories $\vec{z}_{1}(t)$ and $\vec{z}_{2}(t)\equiv\vec{z}_{1}(t)+\vec{\xi}(t)$. Then the largest Lyapunov exponent $\lambda$ can be obtained as (see e.g. Wolf et al. (1985))

[TABLE]

For evaluating this quantity we employed the following numerical procedure Tarkhov et al. (2017): Given a starting point in phase space $\vec{z}_{1}(t_{0})$ we choose a random point $\vec{z}_{2}(t_{0})$ with initial distance $\xi_{0}=\left|\vec{\xi}(0)\right|$ from $\vec{z}_{1}(t_{0})$. We then calculate the time evolution for both trajectories simultaneously by numerical integration. After each time step $l(t)=\log(\xi(t)/\xi_{0})$ is evaluated and recorded, where $\xi(t)=\left|\vec{z}_{1}(t)-\vec{z}_{2}(t)\right|$. When reaching the reset time $T_{r}$ the current Lyapunov exponent is calculated as $\lambda_{t}=l(T_{r})/T_{r}$ and the distance between the trajectories is reset to $\xi_{0}$ by shifting $\vec{z}_{2}(t)$ to $\vec{z}_{2}^{\prime}$ along the vector $\vec{z}_{1}(t)-\vec{z}_{2}(t)$. This procedure is repeated up to a certain time and the largest Lyapunov exponent is calculated as the average over all $\lambda_{t}$ omitting the values up to a time $t_{\mathrm{min}}$. We omit these first values to ensure that we have reached the basin of attraction which is important to get a sensible estimate of the Lyapunov exponent. An alternative method which is inherently linear in the initial separation is evolving the fundamental matrix $\Phi_{t}(\vec{\zeta})$ together with the equations of motion $\dot{\vec{\zeta}}=\vec{F}(\vec{\zeta})$ Skokos (2010):

[TABLE]

Here, $D_{\vec{\zeta}}\vec{F}(\vec{\zeta})$ denotes the Jacobi matrix of the equations of motion. To ensure differentiability, we work in a (real) six dimensional phase space consisting of the real and imaginary parts $\zeta_{i}^{R,I}$, $i\in\{-1,0,1\}$ of the components of $\vec{\zeta}$. At $t=0$ we set $\vec{\zeta}$ according to the desired initial state and $\Phi=\mathbb{1}_{6x6}$. Having calculated the time evolution of the fundamental matrix $\Phi_{t}$, which contains derivatives of the form $\partial\zeta_{i}^{R/I}(t)/\partial\zeta_{j}^{R/I}$ evaluated along the trajectory $\vec{\zeta}(t)$, we obtain the deviation vector at time $t$ as $\vec{\xi}(t)=\Phi_{t}\vec{\xi}(0)$. Choosing again a random initial deviation vector $\vec{\xi}(0)$, this provides an alternative way to calculate the largest Lyapunov exponent via eq. (9). We checked that both methods yield consistent results within numerical errors.

In Fig. 2 we show the Lyapunov exponents calculated as described above for a raster of initial points in the $\left\{\rho_{0},\varTheta_{s}\right\}$-plane corresponding to the same energy used to generate the Poincaré sections in Fig. 1. We find the same qualitative signatures of classical chaos as for the Poincaré sections in Fig. 1: The integrable system ($r=0$) shows no signs of classical chaos. Increasing $r$, we first observe a mixed phase space where it depends on the initial state whether a trajectory has zero or positive Lyapnuov exponent and thus whether it exhibits chaotic behaviour. Finally, for a perturbation $r$ of the same order as the other contributions to the Hamiltonian (2), we see that chaotic regions clearly dominate the phase space. As for the Poincaré sections, the same qualitative features are found for other combinations of the parameters $q$ and $r$.

IV Quantum signatures of chaos

Our aim in this section is to reveal features of the quantum system that can be traced back to chaos in the corresponding classical system and thus may be useful indicators of ergodicity in the quantum system. For this we solve the model (2) by exact diagonalization in Fock space for $N\lesssim 100$, which is feasible as the Hilbert space dimension scales as $N^{2}$.

IV.1 Level statistics

A well-established indicator for chaos – or, more precisely, for randomness – in quantum mechanical systems is the distribution of the energy eigenvalue spacings of the Hamiltonian Mehta (2004).

If a regular system is perturbed, formerly allowed level crossings (e.g. due to symmetries in the unperturbed system) will be avoided. This level repulsion changes the distribution of energy level spacings: In the unperturbed case, as level crossings are allowed, the distribution is peaked at zero and follows a Poisson distribution $P(s)=e^{-s}$. Now with level repulsion, small spacings are suppressed and the distribution of energy level spacings is best approximated by a Wigner distribution $P(s)=\pi/2\ s\ e^{-\pi s^{2}/4}$, which is what is also obtained for random matrices Brody et al. (1981); Graß et al. (2013).

As the Hamiltonian (2) possesses a $\mathbb{Z}_{2}$ symmetry consisting in the invariance under exchange of the state labels $1$ and $-1$, i.e. under $m\mapsto-m$, we first transform to an eigenbasis of this symmetry such that the Hamiltonian matrix becomes block diagonal and calculate the spectrum for each block separately. In addition, we unfold the spectrum in each block to remove the influence of the level density (see e.g. Gubin and F. Santos (2012)) before finally recombining the statistics from the different blocks. The resulting distribution of level spacings $s$ is visualised in Fig. 3 for the values $r=0$ and $r=0.5$ where also the Poisson (dashed green) and Wigner (dotted blue) distributions are shown. To get a quantitative estimate of how strongly chaotic the Hamiltonian is we fit the so-called Brody distribution $P_{b}(s)$ to the distributions (solid red). This function interpolates between the Poisson ($b=0$) and the Wigner distribution ($b=1$) and is given by Brody et al. (1981)

[TABLE]

In Fig. 3 we see that level statistics nicely reflect the classical analysis: For $r=0$ we essentially get a Poisson distribution and for $r=0.5$ the Brody parameter $b\approx 0.6$ suggests that the system is strongly chaotic as expected from the analysis of the classical phase space. The observed $b<1$ is most likely due to the presence of islands of stability at other energies than the one shown in Fig. 1(c). The level spacing distributions shown in Fig. 3 are averaged over the full spectrum. A more thorough analysis would require an energy resolved analysis of the spectral statistics, for which larger Hilbert space dimension would be necessary in order to keep statistical fluctuations small.

A drawback of level statistics as an indicator of quantum chaos is that they only reflect global properties (independent of initial states) and thus do not allow to distinguish the relaxation dynamics of initial states corresponding to different regions of the classical phase space. Also, they are hardly accessible experimentally, which motivates us to turn to experimentally measurable quantities that allow us to resolve the phase space structure, namely Husimi distributions and OTOCs, in the following.

IV.2 Dynamics of Husimi distributions

Husimi distributions allow the visualization of quantum states in a quantum phase space, which is convenient for drawing analogies to the classical picture. In the following, we focus on the case of a classically mixed phase space ($r=0.15$) and consider as initial states coherent states located in the classically chaotic and regular region, respectively, see Fig. 4(a). In order to get an intuition for the differences between these two cases, we study the evolution of the Husimi distribution of each state. The Husimi distribution of a state $\left|\psi\right\rangle$ is a phase space representation and is defined as

[TABLE]

for coherent states $\{\left|\alpha\right\rangle\}$. To obtain the Husimi distribution of a state $\left|\psi\right\rangle$ numerically we would need to calculate the overlap of $\left|\psi\right\rangle$ with each element of a four dimensional grid of coherent states $\{\left|\rho_{0},\varTheta_{s},m,\varTheta_{m}\right\rangle\}$. But as we ultimately want to obtain a two-dimensional phase space representation of $\left|\psi\right\rangle$ that we can compare to the Poincaré sections (see Sec. III), we fix $\varTheta_{m}=0$ and, for each value of $\rho_{0}$ and $\varTheta_{s}$, sum the projections over a grid of $m$-values. The Husimi distributions of two initial coherent states and their time evolution are shown in Fig. 4.

In the limit of large $N$ the size of a coherent state represented in this way would shrink to zero as $1/\sqrt{N}$. Here we use a moderate value of $N=100$ where the extension of the distribution in phase space is still large compared to the classical phase space features we wish to resolve. In fact, the width, i.e. the quantum fluctuations, are so large that both initial states overlap with classically regular as well as chaotic regions of phase space. For the following discussion we should keep in mind that this quantum mechanical smearing over the classical phase space is what fundamentally limits the applicability of semi-classical methods and leads to their breakdown at long times. This point is further discussed in Sec. IV.3 and in the appendix.

Figures 4(b) and (c) show the Husimi distributions of the time-evolved state $\left|\psi(t)\right\rangle$ at time $t=10$. We find that despite the relatively large initial fluctuations qualitative differences between the two cases are clearly visible. The ’regular’ state remains compact in phase space and still resembles a well-localized wave packet while the ’chaotic’ state is distributed across almost the whole classically allowed phase space. This confirms the expectation that the Ehrenfest time, i.e. the time at which quantum interference effects become relevant and the semi-classical wave-packet description breaks down Rammensee et al. (2018); Hummel et al. (2019); Chávez-Carlos et al. (2019); Schubert et al. (2012); Tomsovic et al. (2018); Kurchan (2018); Scaffidi and Altman (2019), is much shorter for chaotic ($\sim\log(N)$) than for regular dynamics ($\sqrt{N}$) Pappalardi et al. (2018).

IV.3 OTOCs

As discussed in Sec. I, semi-classical arguments suggest that OTOCs of the form $C(t)=\langle[\hat{W}(t),\hat{V}(0)]^{\dagger}[\hat{W}(t),\hat{V}(0)]\rangle$ should grow exponentially at short times for systems that exhibit chaotic dynamics in the classical limit. This exponential regime is preceded by an initial powerlaw growth up to the so-called dephasing time and extends up to the Ehrenfest time, or scrambling time. The extension of the window of exponential growth in general depends on the chosen observables and initial states Pappalardi et al. (2019). In the following we discuss to what extent this picture applies to OTOCs with respect to initial coherent states in the three-mode model under study. For models with collective interactions Rozenbaum et al. (2017); Schmitt et al. (2019); Chávez-Carlos et al. (2019); Lewis-Swan et al. (2019); Marino and Rey (2019); Sieberer et al. (2019) the classical limit is well-defined and corresponds to the limit of large particle number. We note that OTOCs have also been studied for systems with local interactions and low-dimensional local degrees of freedom, like spin systems. There, no meaningful classical limit exists and exponential growth is generally not expected Nahum et al. (2018); von Keyserlingk et al. (2018); Xu and Swingle (2020); Khemani et al. (2018). Nevertheless, in this context OTOCs have proven to be useful measures for the spatial spreading of correlations Hosur et al. (2016); Bohrdt et al. (2017); Roberts and Swingle (2016) or diagnosing many-body localization Huang et al. (2017); Swingle and Chowdhury (2017), to give some examples. Here, instead, we focus on a system of collectively interacting bosons and ask whether signatures of the limiting classical dynamics can be observed at finite $N$.

We discuss the time evolution of the squared commutator specifically focusing on the case of $\hat{V}=\hat{W}=\hat{\rho}_{0}\equiv\hat{N}_{0}/N$, and keeping the same parameters and initial coherent states as in the previous section. The evolution of $C(t)$ is shown in Fig. 5(a). The OTOC of the initial state localized in a regular region remains small initially, which is consistent with the regular state evolving quasi-periodically for a long time as observed for the evolution of the Husimi distribution. The OTOC of the chaotic initial state rises extremely fast in comparison. At short times, the onset of the expected exponential behavior is observed in the chaotic case while it is absent in the regular case (inset). At late times, corresponding to the Heisenberg or recurrence time given by the typical inverse level spacing Rammensee et al. (2018); Kurchan (2018), the OTOCs reach a saturation value in both cases.

Next, we study the short-time dynamics (before the Ehrenfest time) in more detail. In this regime we can employ the semi-classical truncated Wigner approximation (TWA), allowing us to simulate arbitrarily large $N$. The squared commutator can be evaluated in TWA through a phase space average over the Poisson bracket Chávez-Carlos et al. (2019); Lewis-Swan et al. (2019); Schuckert and Knap (2019); Pappalardi et al. (2019); Schmitt et al. (2019); Rozenbaum et al. (2017); Cotler et al. (2018) (see the appendix A for details)

[TABLE]

where $\langle\dots\rangle_{W}$ denotes the average over the Wigner function of the initial coherent state. The above expression is based on an expansion in orders of the size of phase space cells $\hbar_{\text{eff}}\sim 1/N$ to order $\hbar_{\text{eff}}^{2}$. We first validate this approach by comparing to exact diagonalization results for $N=100$ in the inset of Fig. 5(a). This shows that while in the regular case, the OTOC is reproduced perfectly up to rather long times, it starts to deviate after the short-time growth regime in the chaotic case. The short-time behavior is always captured well. We note that we had to rescale the TWA results by a factor of order $1$, which depends on the initial state and Hamiltonian parameters but not on the atom number $N$, in order to match the exact diagonalization result. Further benchmark simulations including single-time observables are provided in appendix A showing agreement of TWA with exact diagonalization to much longer times than observed for the OTOCs. We can thus confidently use TWA to explore the short-time dynamics in the large $N$ regime.

As can be seen from Fig. 5(b), particle numbers of $N\gtrsim 10^{4}$ show a more pronounced exponential growth regime at short times. However, this regime does not extend to longer times if we increase the particle number further. Given the observation of exponential growth in the dynamics of the OTOC one may now want to compare the growth rate to the classical Lyapunov exponents shown in Fig. 2. Extracting the growth rate by an exponential fit in Fig. 5(b) we found large deviations from the classical Lyapunov exponent. We emphasize that an agreement between the two is not expected in general. As follows from Eq. (12), the specific OTOC we consider probes the stability in specific phase space directions ($\zeta_{0}^{R}$ and $\zeta_{0}^{I}$) with respect to initial changes in another direction ($\zeta_{0}^{I}$), see Appendix A for details. This corresponds to two specific elements of the fundamental matrix $\Phi$, while the largest Lyapunov exponent is related to the largest eigenvalue of this matrix. As described in Ref. Tarkhov et al. (2017), a derivative such as $\partial\zeta_{0}^{R}(t)/\partial\zeta_{0}^{I}(0)$ evaluated with respect to a single initial condition will in general show erratic behavior including oscillations and at most its envelope is expected to grow exponentially. Exponential growth with the largest Lyapunov exponent $\lambda$ can only be expected after averaging over all initial states of an energy shell. This averaging is what is accomplished by the re-set procedure we use in Sec. III for calculating the Lyapunov exponents. In conclusion, for systems with higher dimensional phase spaces ($d>2$), OTOCs of specific observables and initial coherent states may show clear differences between chaotic and regular initial states but in general do not reveal the largest Lyapunov exponent of the limiting classical dynamics.

We also note that OTOCs of initially commuting operators always show powerlaw growth at short times Pappalardi et al. (2019). The discussions of quantum vs. classical Lyapunov exponents one finds in the literature typically consider cases where the squared commutator is finite at $t=0$ allowing exponential growth to start immediately Prakash and Lakshminarayan (2019); Khemani et al. (2018); Rozenbaum et al. (2017); Schmitt et al. (2019); Stanford (2016).

We calculated $C(t)$ for a large range of other initial conditions and choices of the parameters $q$ and $r$ and evaluated a variety of other OTOCs, for example using spin operators such as $\hat{S}_{z}$ and $\hat{S}_{x}=[\hat{a}_{0}^{\dagger}(\hat{a}_{1}+\hat{a}_{-1})+\rm{h.c.}]/\sqrt{2}$. While the overall features of initial growth and saturation at long times are the same in all cases, we cannot always clearly identify exponential growth at short times for classically chaotic initial states which we attribute to the fact that we probe specific elements of the stability matrix, which may show arbitrary behavior in general Tarkhov et al. (2017), and that the time window between the initial powerlaw growth and the Ehrenfest time may be small. For small values of $r$ we consistently found that if the initial state overlaps with the position of the separatrix of the (regular) classical phase space, the OTOC $C(t)$ grows much faster than if it lies close to the stable fixed points. We conclude that despite OTOC growth not being directly related to the largest Lyapunov exponent of the corresponding classical dynamics, OTOCs show clear differences between initial states localized in classically chaotic and regular regions of phase space, thus qualifying as indicators of quantum chaos.

In order to simulate dynamics past the Ehrenfest time, and in particular to study the asymptotic behavior of OTOCs in the long-time regime we need to resort to exact diagonalization. An exciting alternative is to measure OTOCs in quantum simulation experiments with spinor condensates, discussed in the following, where particle numbers up to $10^{4}$ are reached easily.

V Experimental realization

In this section we discuss a protocol for measuring OTOCs experimentally based on time reversing the system dynamics and propose an implementation using a spinor BEC of rubidium atoms.

The most intuitive protocol for testing reversibility and chaos in a quantum system is the following: Prepare an initial state, evolve it to time $t$, apply a perturbation $\hat{W}$, evolve for another period $t$ under the sign-reversed Hamiltonian, and then measure an observable $\hat{V}$. This protocol results in the measurement of $\langle\hat{W}(t)^{\dagger}\hat{V}\hat{W}(t)\rangle$. If the initial state is an eigenstate of $\hat{V}$, then the measured quantity has the desired form of the out-of-time-order part of the squared commutator Gärttner et al. (2017). Note, however, that in this protocol $\hat{W}$ is the unitary operator that induces the applied perturbation. In order to make the connection to Lyapunov exponents we need $\hat{W}$ to be a Hermitian operator. For spin-$1/2$ chains one can for example choose $\hat{W}$ to be a spin flip on one of the spins, $\hat{\sigma}_{i}^{x}$, which is Hermitian (and unitary). In the case of bosonic systems, however, this simplification cannot be made.

A possible ansatz for solving this problem is to consider small perturbations and to extract the resulting signal to low order in the perturbation strength. Mathematically, we consider $\hat{W}(0)=\exp[-i\phi\hat{A}]$, with a Hermitian operator $\hat{A}$. To second order in $\phi$ one obtains

[TABLE]

where the contribution quadratic in $\phi$ contains an OTOC of Hermitian operators: $\hat{A}$ acts at time $t$, then $\hat{V}$ acts at time [math] followed by $\hat{A}$ acting at time $t$. This ordering is not compatible with the ordering of operators on a Schwinger-Keldysh contour where the ordering is from small to long and back to small times, hence the name OTOC. In general, it is not straight forward to separate the OTOC part from the additional terms that appear. However, in the case where the initial state $\left|\psi_{0}\right\rangle$ is an eigenstate of the operator $\hat{V}$ matters simplify considerably and the squared commutator can be obtained from measurements of the quadratic response after time reversal as we discuss in the following.

For $\hat{V}\left|\psi_{0}\right\rangle=\Lambda\left|\psi_{0}\right\rangle$ the squared commutator becomes (expectation values are always with respect to $\left|\psi_{0}\right\rangle$)

[TABLE]

For the response of the n$th$ moment of $\hat{V}$ the commutator term [linear term in $\phi$ in Eq. (13)] vanishes such that up to second order in $\phi$ one obtains

[TABLE]

where we have defined the quadratic response $\Gamma^{(2)}_{\hat{O}}$ of an operator $\hat{O}$. Thus the squared commutator can be obtained by measuring the quadratic response of the first and second moment of $\hat{V}$ in the time-reversal protocol described above as

[TABLE]

In the case at hand one could for example initialize the system in the state where all particles are in the $m_{F}=0$ mode and measure $\hat{V}=\hat{N}_{0}$.

We now turn to a concrete experimental implementation of the model studied above. As mentioned in Sec. I, it can be straightforwardly realized using a BEC of 87Rb atoms in the $F=1$ hyperfine manifold in a tightly confining trap. The integrable case $r=0$ has been realized experimentally and its dynamics studied by a number of groups Stamper-Kurn and Ueda (2013); Gerving et al. (2012); Hamley et al. (2012); Lange et al. (2018); Kunkel et al. (2018); Linnemann et al. (2016). The integrability breaking term can be implemented by continuously applying an rf driving field where $r$ corresponds to the Rabi frequency induced by the drive. rf-driving is commonly used to implement spin rotations that are fast compared to the interacting dynamics in order to read out different components of the spin Kunkel et al. (2019). By reducing the modulation amplitude of the magnetic field the Rabi frequency can be adjusted such that $r$ is widely tunable. The interaction term in the Hamiltonian results from s-wave collisions and is typically on the order of $gN/2\pi\approx-2\,$Hz. The parameter $q$ can be adjusted by microwave dressing to be of the same order Kunkel et al. (2018).

The challenge for experiments lies in measuring the OTOC, which involves time-reversing the unitary dynamics of the system as illustrated by the protocols discussed above. Indeed, for 87Rb the nature of the collisional interactions between the different hyperfine components allows one to do exactly this. In the $F=2$ hyperfine manifold the sign of $g$ in the spin-changing collision term is opposite to $F=1$ and significantly larger. The other parameters of the Hamiltonian can be adjusted freely through external control fields. Thus, reversal of the dynamics can be accomplished by transferring the BEC into the $F=2$ manifold and running the dynamics under the sign-reversed Hamiltonian for a time that has to be adjusted in order to compensate for the different magnitude of the Hamiltonian parameters. The additional two levels ($m_{F}=\pm 2$) in the $F=2$ manifold can be tuned out of resonance by adjusting the magnetic field and thus the quadratic Zeeman shift, or by off-resonant microwave dressing, resulting in effective three-mode dynamics.

Finally, we discuss a number of recent proposals for measuring OTOCs. We tried to adapt them to the proposed implementation, however, all of them either measure OTOCs with respect to thermal (or even infinite temperature) states or have technical requirements that seem extremely challenging for spinor BECs. Protocols using randomized measurements Vermersch et al. (2019) only allow to determine OTOCs of the maximally mixed state. The scheme proposed in Ref. Bohrdt et al. (2017) requires preparing two identical copies of the system and detecting the atom number parity, while Yao et al. (2016) additionally requires an ancillary systems that controls the sign of the Hamiltonian and similarly for Zhu et al. (2016) which is designed to detect OTOCs of thermal states.

VI Conclusions and outlook

This study shows that a three-mode bosonic system with spin-changing collisions and tunable integrability breaking rf-coupling term can serve as a showcase system for studying classical and quantum chaos. In the classical limit the system transitions from regular to chaotic phase space dynamics as the strength of the rf-coupling is increased from zero to $\sim 1$ in units of the characteristic energy scale of the interactions. This transition is also confirmed by the Lyapunov exponents which are positive in the emerging chaotic regions of phase space. The statistics of eigenenergies, a traditional indicator of quantum chaos, behaves as expected and changes from Poissonian to Wigner-Dyson as the integrability breaking term is tuned. For OTOCs, which have been proposed as faithful indicators for quantum chaos in the semi-classical regime, we find faster growth for initial states close to the separatrix where the classical dynamics becomes chaotic even for small integrability breaking. The characteristic exponential growth that is expected due to the connection between the OTOC and the Lyapunov exponent in the classical limit is visible in an appreciable time window for large particle numbers ($\gtrsim 10^{3}$). A direct comparison of the resulting growth exponent to the largest Lyapunov exponent defined in the standard way is not meaningful here, which can be understood by considering the semi-classical truncated Wigner expression (12) for the OTOC. This form shows that an OTOC of specific operators will depend on specific elements of the stability matrix, or fundamental matrix, of the classical dynamics, which are in general not expected to grow exponentially with the largest Lyapunov exponent.

The considered model can be implemented experimentally, potentially giving access to OTOCs beyond the Heisenberg time for much larger particle numbers than what is reachable by exact diagonalization. Time-reversal may be achieved in BECs of 87Rb atoms by exploiting the collisional properties in different hyperfine states. This enables the measurement of OTOCs with respect to initial states that are eigenstates of the measured operator. Specifically the squared commutator, which has a particularly direct semi-classical meaning, can be extracted for initial coherent states which are easy to prepare experimentally. The protocol for extracting the squared commutator requires measuring the quadratic response with respect to a perturbing generalized rotation. More straightforward protocols for testing the reversibility and sensitivity of the many-body dynamics with respect to perturbations also allow measuring quantities that are of interest for semi-classical physics and quantum chaos Schmitt et al. (2019). Therefore, in the future, we will seek to clarify the meaning of the resulting more general types of OTOCs in the context of operator spreading and thermalization. Intuitively, these objects should still be indicators of ergodicity and emerging operator complexity, and are in addition related to quantum Fisher information Gärttner et al. (2018), enabling connections to entanglement.

Acknowledgements.

We thank Thomas Gasenzer, Benjamin Geiger, Quirin Hummel, Daniel Linnemann, and Markus Oberthaler for discussions. Supported by the DFG Collaborative Research Center SFB1225 (ISOQUANT).

Appendix A Truncated Wigner Approximation

The TWA relies on a phase space representation of quantum states where expectation values of quantum operators can be calculated by evaluating their corresponding Weyl symbols averaged over phase space weighted with the Wigner function where the latter is the Weyl symbol of the density matrix of the initial state. As we are interested in calculating OTOCs, we first note that the Weyl symbol of the commutator $\left[\hat{A},\hat{B}\right]$ is given by

[TABLE]

where $\left\{\dots\right\}_{MB}$ denotes the Moyal bracket, $\left\{\dots\right\}$ the Poisson bracket and $O_{W}$ is the Weyl symbol of the quantum operator $\hat{O}$ Polkovnikov (2010). $\hbar_{\text{eff}}$ quantifies the minimal phase space volume that a state can be confined to due to Heisenberg’s uncertainty relation, which corresponds to the phase space volume of a coherent state, $\hbar_{\text{eff}}\sim 1/N$ in our case. Working directly in the $\zeta$ representation (see beginning of section II) where the real and imaginary parts of $\zeta_{i}$ act as pairs of conjugate variables, the Poisson bracket is defined as

[TABLE]

i.e. we have for $A=\rho_{0}(t)=\left(\zeta_{0}^{R}(t)\right)^{2}+\left(\zeta_{0}^{I}(t)\right)^{2}$ and $B=\rho_{0}(0)$ up to order $\hbar_{\text{eff}}^{2}$

[TABLE]

where $\langle\dots\rangle_{W}$ denotes the phase space average weighted by the Wigner function. In the second step we used that we can always choose $\zeta_{0}$ to be real initially by choosing an appropriate global phase. The derivative ${\partial\rho_{0}(t)}/{\partial\zeta_{0}^{I}}$ is evaluated numerically by evolving two trajectories initially separated by ideally infinitely small $d_{0}$ in the imaginary part of the $\zeta_{0}$ component and evaluating the difference in $\rho_{0}$ at each time step numerically integrating the mean field e.o.m. (4). Alternatively, we can rewrite

[TABLE]

where now the derivatives can be evaluated as elements of the fundamental matrix as described in Sec. III.2. We again checked that the two methods yield the same results for small enough initial deviations $d_{0}$ in the case of the two trajectory calculation (or equivalently, for short enough times). The Wigner function of the initial coherent state is well approximated by a Gaussian with mean and variance matching the expectation value and variance of the initial state w.r.t. the quantum operators. Technically, we sample from $\left|(x_{1}+ix_{2})/2,N,(x_{3}+ix_{4})/2\right\rangle$ where all $x_{i}$ are independently drawn from a Gaussian distribution with zero mean and standard deviation $1/\sqrt{N}$ to approximate the Wigner function of the initial state $\left|0,N,0\right\rangle$. To get the Wigner function of an arbitrary initial coherent state, we rotate the sampled states using spin operators as generators of the rotation.

We checked that the sampling of arbitrary initial states as described above works by comparing the expectation values and variances of different spin operators calculated within TWA to the exact quantum mechanical calculations for up to $N=100$ particles, see Fig. 6. We also checked that in the case of the two trajectory calculation $d_{0}$ is small enough s.t. the difference quotient is independent of $d_{0}$ and matches the quantum calculations after appropriately scaling the amplitude (see inset of Fig. 5). We can confirm that the TWA results match the exact quantum calculations up to $t\approx 2.5$ for $N=100$ particles and expect that the approximation holds up to longer times with increasing $N$ as $t\sim\log N$.