Quantum chaos in a Bose-Hubbard dimer with modulated tunnelling

TL;DR

This paper investigates quantum chaos in a Bose-Hubbard dimer with modulated tunnelling, using numerical simulations and statistical measures to distinguish chaotic from regular quantum dynamics, and demonstrates the effectiveness of the truncated Wigner approximation.

Contribution

It introduces the statistical distance between number distributions as a new measure to identify chaos in quantum dynamics of the Bose-Hubbard dimer.

Findings

Statistical distance effectively distinguishes chaotic from regular quantum states.

Truncated Wigner approximation accurately captures beyond mean-field dynamics.

Chaotic regions exhibit irregular Wigner functions, yet are well-represented by the method.

Abstract

In the large-$N$, classical limit, the Bose-Hubbard dimer undergoes a transition to chaos when its tunnelling rate is modulated in time. We use exact and approximate numerical simulations to determine the features of the dynamically evolving state that are correlated with the presence of chaos in the classical limit. We propose the statistical distance between initially similar number distributions as a reliable measure to distinguish regular from chaotic behaviour in the quantum dynamics. Besides being experimentally accessible, number distributions can be efficiently reconstructed numerically from binned phase-space trajectories in a truncated Wigner approximation. Although the evolving Wigner function becomes very irregular in the chaotic regions, the truncated Wigner method is nevertheless able to capture accurately the beyond mean-field dynamics.

Figures

Equations

$$\hat{H}=\hslash U\left({\hat{b}}_1^{†}{\hat{b}}_1^{†}{\hat{b}}_1{\hat{b}}_1+{\hat{b}}_2^{†}{\hat{b}}_2^{†}{\hat{b}}_2{\hat{b}}_2\right)-\hslash J\left({\hat{b}}_1^{†}{\hat{b}}_2+{\hat{b}}_2^{†}{\hat{b}}_1\right),$$

$$x$$

$$y$$

$$z$$

$$\frac{\mathrm{d}z}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\varphi }{\mathrm{d}t}$$

$$\lambda =\underset{t\to \infty }{\lim }\underset{{\mathbf{z}}_0^{\prime }\to {\mathbf{z}}_0}{\lim }\frac{1}{t}\ln \left(\frac{\left|\left|\mathbf{z}(t)-{\mathbf{z}}^{\prime }(t)\right|\right|}{\left|\left|\mathbf{z}(0)-{\mathbf{z}}^{\prime }(0)\right|\right|}\right),$$

$${\hat{H}}_{\text{eff}}$$

$$\displaystyle\quad-\frac{\hbar U}{4}\left[1-\mathcal{J}_{0}\left(\frac{4\mu}{\omega}\right)\right]\Big{(}\hat{b}_{1}^{\dagger}\hat{b}_{1}^{\dagger}\hat{b}_{2}\hat{b}_{2}+\hat{b}_{2}^{\dagger}\hat{b}_{2}^{\dagger}\hat{b}_{1}\hat{b}_{1}$$

$$\displaystyle\qquad\qquad-4\hat{b}_{1}^{\dagger}\hat{b}_{1}\hat{b}_{2}^{\dagger}\hat{b}_{2}\Big{)}-\hbar J_{0}\left(\hat{b}_{1}^{\dagger}\hat{b}_{2}+\hat{b}_{2}^{\dagger}\hat{b}_{1}\right)$$

$${\hat{\rho }}_{\text{red}}=\left(\begin{array}{cc}⟨{\hat{b}}_1^{†}{\hat{b}}_1⟩& ⟨{\hat{b}}_1^{†}{\hat{b}}_2⟩\\ ⟨{\hat{b}}_2^{†}{\hat{b}}_1⟩& ⟨{\hat{b}}_2^{†}{\hat{b}}_2⟩\end{array}\right).$$

$${\lambda }_{\text{max}}$$

$$|⟨\hat{\mathbf{J}}⟩|$$

$$\frac{\mathrm{d}{\beta }_1}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\beta }_2}{\mathrm{d}t}$$

$$\frac{1}{2}⟨{\hat{b}}_i^{†}{\hat{b}}_j+{\hat{b}}_j{\hat{b}}_i^{†}⟩=\int {\beta }_i^{\ast }{\beta }_{j}W({\vec{\beta }}^{\ast },\vec{\beta }){\mathrm{d}}^4\beta .$$

$$O(\psi ,{\psi }^{\prime })=\left|⟨\psi |{\psi }^{\prime }⟩\right|=\left|\sum _{n,m}{c}_n^{\ast }{d}_m⟨n|m⟩\right|=\left|\sum _n{c}_n^{\ast }{d}_n\right|,$$

$$B(\psi ,{\psi }^{\prime })=\sum _n\sqrt{P_n\widetilde{P_n}}=\sum _n\left|c_n^{\ast }{d}_n\right|.$$

$$D_B(\psi ,{\psi }^{\prime })=-\ln \left[B(\psi ,{\psi }^{\prime })\right],$$

$$n\le {\left|{\beta }_1\right|}^2<n+1.$$

$$\frac{{\sigma }_p}{\bar{p}}$$

$${\hat{J}}_x$$

$${\hat{J}}_z$$

$${\hat{H}}_0$$

$${\hat{H}}_1$$

$${\hat{H}}_{\text{eff}}=\frac{1}{T}{\int }_0^T{e}^{iA(t){\hat{J}}_x}{\hat{H}}_0{e}^{-iA(t){\hat{J}}_x}\mathrm{d}t,$$

$$e^{iA(t){\hat{J}}_x}$$

$$=\frac{\cos (2A(t))}{2}\left({\hat{J}}_z^2-{\hat{J}}_y^2\right)$$

$$+\frac{\sin (2A(t))}{2}\left({\hat{J}}_y{\hat{J}}_z+{\hat{J}}_z{\hat{J}}_y\right)$$

$$+\frac{1}{2}\left({\hat{J}}_z^2+{\hat{J}}_y^2\right)$$

$${\int }_0^T\cos (2A(t))\mathrm{d}t$$

$${\int }_0^T\sin (2A(t))\mathrm{d}t$$

Full text

quantum chaos; ultracold atoms

I Introduction

Chaos is characterised by an exponential sensitivity to small perturbations in initial conditions (hyperbolicity), non-integrable dynamics, and a widely-traversed phase space (ergodicity) Wimberger (2014). Despite its ubiquity in the classical macroscopic world, chaos is not intrinsic to the underlying quantum description, except where an effective nonlinearity is introduced through, e.g., taking a classical limit or mean-field approximation, or conditioning the dynamics upon the results of a measurement Eastman et al. (2017). Nevertheless, the behaviour of quantum states in classically chaotic regimes, referred to as ‘quantum chaos’, can mimic features of classical chaos, at least for short times. Quantum chaos is thus of intrinsic interest for the insight it brings to the quantum-classical crossover Haake (2010); Gutzwiller (1992), in addition to its connections with ergodicity and thermalisation Altland and Haake (2012).

The sensitivity to initial conditions in classical chaos is quantified by Lyapunov exponents: a positive Lyapunov exponent indicates an exponential growth in the separation in phase space between trajectories that are initially close together. A direct implementation of this measure in quantum systems is problematic due to the intrinsic quantum uncertainty: it is not possible deterministically to specify arbitrarily precise initial conditions in phase space. A naive alternative would be the overlap of initially similar quantum states; however this overlap is time-invariant under unitary evolution. For this reason, the overlap of initially identical distributions but under slightly perturbed Hamiltonians has been used to distinguish between quantum chaos and regular dynamics Haake (2010). States with large support on regular regions of phase space will retain close to unitary overlap under a small Hamiltonian perturbation, whereas states that are well-confined to chaotic regions will exhibit an exponential decay of overlap to around order $1/N$ Haake (2010), where $N$ is the number of particles. This decay may be used to extract effective Lyapunov exponents Cucchietti et al. (2002).

However, quantum-state overlap is not directly accessible experimentally. Even theoretically, the quantum state becomes intractable to calculate or even specify exactly for moderately sized systems, due to the exponential scaling of Hilbert-space dimension on the number of modes.

In this paper, we propose and benchmark the use of the overlap of number distributions, rather than the quantum state itself, as a measure of quantum chaos. Number distributions are experimentally accessible in ultra-cold atom experiments, and can be calculated without knowing the full quantum state, such as through the binning of trajectories in a truncated Wigner simulation Lewis-Swan et al. (2016). In particular, we calculate the Bhattacharyya distance Bhattacharyya (1943) between number distributions in a driven, double-well Bose-Einstein condensate, bench-marking approximate approaches against exact calculations. We find that the Bhattacharyya distance – even when calculated approximately using binned truncated Wigner trajectories – is a reliable indicator of quantum chaos that is consistent with other measures. The advantage of this method is that it scales well with system size and so may be applied to a wide range of systems for which exact calculations are intractable. Our results also provide further insight into the nature and scope of quantum chaos in a simple model that can be implemented and probed experimentally Tomkovič et al. (2017).

We note that an alternative solution is possible for continuously monitored quantum systems Eastman et al. (2017); Pokharel et al. (2018); Eastman et al. (2019), where individual stochastic quantum trajectories, representing the dynamics conditioned on measurement, can be used to derive effective quantum Lyapunov exponents. In this work, we study exclusively the isolated dynamics of the two-site lattice, where such an approach based on continuous measurement readout is not available.

Distributions over phase space provide a powerful way to connect the evolving quantum state to the dynamical structure of the corresponding classical system. For example, the Husimi $Q$-function, besides providing a visualisation of the state Mahmud et al. (2005), is the basis for an approximate treatment of the quantum dynamics through an ensemble of classical trajectories Trimborn et al. (2009); Hennig et al. (2012). Our use of a Wigner function aims to improve on this approach, since to obtain classical Liouvillian evolution, the truncation of terms in the evolution equation is less severe in the Wigner case than in the $Q$-function case. Although the truncated Wigner method has been used extensively in ultracold atoms, it has not, to our knowledge, been used to study chaotic dynamics in quantum systems.

We begin, in section II by introducing the double-well system and discussing the chaotic dynamics that arises in a semiclassical analysis when the tunnelling rate is modulated in time. In section III we consider the quantum dynamics, calculated exactly and with the truncated Wigner method. Finally in section IV we introduce the Bhattacharyya distance as a measure of chaos in quantum systems and test its performance when used with the truncated Wigner representation.

II Modulated Bose-Hubbard dimer

For a sufficiently tight trapping potential, a Bose-Einstein condensate in a double-well potential Milburn et al. (1997a) may be modelled as a two-site Bose-Hubbard model, or quantum dimer:

[TABLE]

where $J$ is the tunnel-coupling rate and $U$ is the strength of the on-site interaction term, which arises from interparticle collisions treated in the low-temperature $s$-wave limit. The canonical boson operators obey the commutation relations $[\hat{b}_{j},\hat{b}_{k}^{\dagger}]=\delta_{jk}$. To introduce chaos into this system, we will vary the tunnelling rate in time: $J(t)=J_{0}+\mu\cos{(\omega t)}$, which could be achieved, for example, through modulation of the barrier height Milburn et al. (1997b).

Even for thousands of particles, the relevant Hilbert space of the system is small enough to allow for direct calculation of the dynamical wavefunction, affording the opportunity to benchmark approximate approaches against exact solutions; yet the interplay between coherent tunnelling and particle interactions leads to rich many-body physics. This physics is reflected in the nonlinear behaviour of the corresponding semiclassical model, which can undergo macroscopic self-trapping and, for appropriate modulation of parameters, chaotic dynamics.

The semiclassical, or mean-field, model describes the dynamics of the order parameter under the assumption that two-body correlations factorise, $\langle AB\rangle=\langle A\rangle\langle B\rangle$. It is justifiable in the large-particle limit where fluctuations around the mean value of observables are small. It is convenient to define ‘Bloch-sphere’ observables:

[TABLE]

which, in the semiclassical limit, satisfy $x^{2}+y^{2}+z^{2}=N^{2}/4$. The dynamics on the Bloch sphere can be mapped to a plane $(z,\phi)$, where $z$ is the population difference between the two sites and $\phi=-\arg{(x+iy)}$ is the relative phase. The semiclassical dynamics is then given by the coupled equations:

[TABLE]

Without modulation of the parameters, this system is integrable, with solutions given by Jacobi elliptic functions Milburn et al. (1997a). The nonlinear dynamics is governed by the ratio of interaction energy to the tunnelling rate, $C=UN/J_{0}$. As $C$ is increased through the critical value $C=1$, the stable centre at $(z,\phi)=(0,0)$ undergoes a pitchfork bifurcation. The separatrix that emerges from the resultant saddle point encloses each of the two new stable centres and prevents trajectories in one enclosed region from traversing to the other. This nonlinear suppression of tunnelling, known as macroscopic self-trapping, has been observed experimentally Albiez et al. (2005).

Modulating the tunnelling rate breaks integrability and gives rise to the possibility of chaos. For example, at the critical value $C=1$, a chaotic sea emerges in the region of the bifurcating point (see Fig. 1) and continues to persist at $C>1$ around the separatrix Milburn et al. (1997b); Mourik et al. (2018). We note that integrability may be broken in other ways, such as modulation of the interaction strength Gong et al. (2009); Watanabe and Mäkelä (2012) or the relative height of the wells Morales-Molina and Gong (2008).

Chaos in phase space can be illustrated by stroboscopic Poincaré sections, which are constructed by evolving and periodically plotting a selection of semiclassical trajectories at intervals of the driving period Haake (2010). Poincaré sections at the bifurcation point are shown in Fig. 1a and Fig. 1b, without and with driving, respectively. In Fig. 1b a chaotic sea emerges around the bifurcating fixed point, within which sit several Kolmogorov-Arnold-Moser (KAM) islands Haake (2010).

The maximal Lyapunov exponent at phase-space coordinate $\bm{z}(0)$ is defined as

[TABLE]

where the initial values of the two phase-space trajectories $\bm{z}(0)$ and $\bm{z}^{\prime}(0)$ should in principle be made vanishingly close to each other. In practice, we implement this numerically by choosing a small but finite difference and simulate for a time that is long enough for the maximal exponent to dominate the rate of divergence (but small enough that this divergence is not limited by phase-space boundaries) Chianca and Olsen (2011); Wimberger (2014). Fig. 1c was generated by evolving pairs of perturbed points distributed across phase space and performing a linear fit on their logarithmic relative divergence, $\ln{\left(\frac{\left|\left|\bm{z}(t)-\bm{z}^{\prime}(t)\right|\right|}{\left|\left|\bm{z}_{0}-\bm{z}^{\prime}_{0}\right|\right|}\right)}$ against time, $t$. The initial perturbation magnitude was $\left|z-z^{\prime}\right|=10^{-4}$ and the trajectories were evolved for 20 driving periods. The evolution time of the finite-time Lyapunov exponents was small enough to avoid saturation of the logarithmic relative divergence. The darker regions of Fig. 1c, which represent large positive Lyapunov exponents, correspond well with the chaotic regions of phase space shown in the Poincaré section, Fig. 1b. The white regions in Fig. 1c correspond to non-positive Lyapunov exponents.

The extent of chaos in the modulated dimer can be controlled via the driving amplitude, $\mu$, and frequency, $\omega$. Fig. 2 quantifies the fraction of phase space occupied by chaotic regions. For this plot, we calculated the finite-time Lyapunov exponents for 1600 trajectories distributed across the phase space and classified each phase-space point as chaotic if its Lyapunov exponent was greater than the largest Lyapunov exponent within the unmodulated phase space. The percentage of phase space that was chaotic was then determined and plotted for varying modulation amplitude and frequency. We note that this approach does not produce exact Lyapunov exponents, yet is sufficient to distinguish regular and chaotic trajectories to a high degree of accuracy in all tested parameter regimes.

Even though the modulation is a requirement for chaotic behaviour in this system, driving at high frequency tends to diminish the extent of the chaos. We illustrate this Fig. 2 by plotting several representative Poincaré sections highlighting the different regimes. For example, a sufficiently rapid modulation ($\omega\gtrsim 6J_{0}$) will suppress chaos entirely, whereas for low modulation frequency ($\omega\sim 0.5J_{0}$), widespread chaos persists for modulation amplitudes as high as $\mu=100J_{0}$.

Consistent with the disappearance of chaos at high modulation frequency observed in Fig. 2, an effective static integrable Hamiltonian Kohler et al. (2005); Gong et al. (2009) can be derived under the assumption that the modulation frequency is much larger than other frequency scales in the system:

[TABLE]

where $\mathcal{J}_{0}$ is the zeroth order Bessel function of the first kind (calculations in App. A). As shown in Figs. 3 and 4, the Poincaré sections corresponding to the time-periodic and effective Hamiltonians become very similar for sufficiently large modulation frequency ($\omega=10$).

III Quantum dynamics

For the two-mode system, the semiclassical results above can be compared to the exact quantum dynamics, which can be calculated through a Floquet analysis. Such an approach was used, for example, to investigate dynamical tunnelling between regular islands in the semiclassical Poincare section Salmond et al. (2002); Watanabe and Mäkelä (2012).

To make a connection with the classical Poincare sections, we represent the evolving quantum state in phase space through use of the atomic $Q$ function Arecchi et al. (1972); Lee and Feit (1993). Fig. 5 shows the $Q$ function of what are initially Bloch coherent states 111Bloch coherent states Arecchi et al. (1972) are minimal uncertainty states of the dimer under the restriction of fixed total particle number. They may be created experimentally by, for example, forming a condensate of known number in one well, followed by linear tunnelling Tomkovič et al. (2017). centred at three representative points in phase space.

After 20 driving periods, the $Q$ function initially centred in the chaotic region has become very irregular and is smeared out over nearly all of the chaotic sea, reflecting the ergodic nature of the corresponding classical trajectories Mahmud et al. (2005). By contrast, the $Q$ function of the ‘regular 1’ state remains localised on the KAM island on which it was initially placed. The ‘regular 2’ state is smeared out, but in a qualitatively very different way from that seen in the chaotic region; the dependence of the tunnelling oscillation period upon its amplitude causes the distribution to be sheared along the regular classical trajectories (a two-mode analogue of the optical Kerr effect).

When the $Q$ function is no longer well-localised, it is clear that the semiclassical predictions – which neglect fluctuations about the mean – will become unreliable. This neglect of fluctuations in the chaotic regime is rapidly manifest in discrepancies in the predictions for observables, as illustrated in Fig. 6. In the chaotic regime, the semiclassical simulation shows persistent, albeit irregular, tunnelling oscillations in the number difference $\langle\hat{z}\rangle$ whereas in the quantum dynamics this oscillation rapidly decays. The decay in the mean of $z$ is associated with a rise in the variance, which supports the interpretation of the chaotic quantum state as an ensemble of irregular classical trajectories that get out of phase over time. On a much longer time scale, the tunnelling oscillations of the ‘regular 2’ state also collapse, but this is due to the gradual dephasing of regular oscillations at slightly different frequencies.

The fact that a fairly rapid growth of variance occurs in the ‘regular 2’ as well as the chaotic regions indicates that the variance itself is not always a reliable probe of the presence of chaos Tomkovič et al. (2017). Nevertheless, when growth in variance is driven by chaotic behaviour, it does reveal something about the extent of the chaotic region.

For the two-mode system, the spreading of the initially coherent $Q$ function over the Bloch sphere is directly linked with the fragmentation of the condensate Trimborn et al. (2009). Fragmentation occurs when not all particles occupy the same single particle state – as required for a pure BEC – and is quantified by the eigenvalues of the one-particle reduced density matrix Leggett (2001),

[TABLE]

A fragmented condensate is revealed when $\hat{\rho}_{\text{red}}$ has more than one nonzero eigenvalue, i.e., when $\lambda_{\text{max}}<N$.

For the two-mode system, the maximum eigenvalue of $\hat{\rho}_{\text{red}}$ can be written in terms of the length of the Bloch spin vector $|\langle\hat{\bf{J}}\rangle|$ Trimborn et al. (2009), which in turn can be expressed in terms of the variance of the Bloch variables:

[TABLE]

Since fragmentation is linked to the growth in the variance of any of the Bloch variables, it provides a more global picture of the spread of the state over the Bloch sphere than the variance in any one particular variable.

The largest eigenvalue is plotted as a function time in Fig. 7 for the three representative initial conditions, both with and without modulation. The first thing to note is that loss of coherent tunnelling oscillations does not mean total fragmentation: the maximum eigenvalue for the state in the chaotic sea falls only to about $0.9N$, far above the lower limit of $0.5N$. Second, fragmentation is not just a consequence of chaotic behaviour. Indeed, the largest fragmentation occurs in the regular region where the classical dynamics exhibit large-amplitude tunnelling oscillations Hennig et al. (2012) and for which the corresponding $Q$ functions are sheared across a large solid angle on the Bloch sphere. The $Q$ function of the state in the chaotic region – although it spreads out to subtend a similarly large angle on the Bloch sphere – actually corresponds to less fragmentation, since its $Q$ function fills in the intermediate region of phase space and thus maintains a higher level of total average spin. For larger relative nonlinearity $C$ or larger modulation amplitude $\mu$, where the chaotic sea occupies a greater fraction of phase space Trimborn et al. (2009) (see Fig. 2), we see increased fragmentation from the chaotic dynamics.

Third, perhaps the most striking consequence of the modulation is the near total suppression of fragmentation in the $Q$ function of the ‘regular 1’ state. This suppression of fragmentation is due to the localisation enforced by creation of the KAM islands, and is particularly clear in Fig. 7b, in which the regions of high purity (low fragmentation) match closely the centres of the KAM islands seen in in Fig. 1.

Overall then, in regions of phase space dominated by chaotic behaviour, the structures in the Poincare section determine the level of fragmentation in the quantum dynamics. Conversely, if the system is known to be chaotic, then fragmentation in the quantum system could be used to determine the presence of KAM islands or to probe the size of the chaotic sea, even though fragmentation by itself is not an indicator of the presence of chaos.

By choosing system parameters to precisely alter the size of the chaotic sea, it is possible to generate states with arbitrary condensate fraction from initial coherent states on the timescale of few modulation periods.

III.1 Truncated Wigner method

Thus far we have calculated the quantum dynamics through exact calculation of the state vector, with a subsequent mapping to phase space via the $Q$ function. We now use these results to benchmark the truncated Wigner method as a tool for investigating quantum chaos. Not only is the truncated Wigner method an efficient way to simulate many-mode quantum systems (the exact state-vector approach quickly becomes intractable even for just a few modes), it gives a direct access to dynamics in phase space, in which much of the classical analysis of chaotic systems is framed.

In the truncated Wigner method Graham (1973), a Fokker-Planck equation for an approximate Wigner function is obtained by truncation of the full partial differential equation for the Wigner function at second-order. This approximation is justified by the relatively small size of the coefficients of the third-order terms for a system with many particles. The Fokker-Planck equation is then mapped through standard techniques Gardiner and Collett (1985) onto a set of ordinary differential equations for stochastic complex amplitudes, $\beta_{j}$. Although the Wigner function is not positive-definite for all states, a positive approximation to common initial states Olsen and Bradley (2009) allows the low-order moments to be accurately sampled. Any such approximation to the initial state together with the truncation of the third-order terms neglects any negativity of the Wigner function, but gives accurate operator moments up to the collapse timescale of quantum optics Corney and Olsen (2015).

The truncated Wigner equations for the Bose-Hubbard dimer are:

[TABLE]

which are solved for an ensemble of trajectories. Expectation values of symmetrically ordered products can be calculated directly as stochastic averages, for example,

[TABLE]

The truncated Wigner simulations used in this work were implemented with the XMDS2 Dennis et al. (2013) software package.

The truncated Wigner calculation of the expectation value of the two-mode number difference, $z$, is included in Fig. 6. The results for the two states in regular regions are indistinguishable from the exact calculations over this time scale. Some discrepancy is seen for the state in the chaotic region, but only well after the collapse of the tunnelling oscillations.

Binning trajectories allows us to build up a picture of the phase-space distribution without having to calculate the full density matrix. For large $N$, the exact Wigner function on the Bloch sphere is particularly difficult to calculate, as it requires generation of high-order spherical harmonics Trimborn et al. (2009). Fig. 8 shows the binned distribution after 20 driving periods for the three representative states, showing a close correspondence to the $Q$ functions calculated in Fig. 5.

IV Bhattacharyya distance

Sensitivity to initial conditions is the hallmark of chaos in classical systems. Although a direct measure of the divergence of particular phase-space trajectories is not accessible in unitary quantum dynamics, we have seen above that the divergence is reflected in the way that $Q$ and Wigner functions within the chaotic region spread out to cover the whole of the chaotic sea. Moreover, they evolve in a very irregular way that contrasts markedly with the regular shearing, for example, that occurs for the ‘regular 2’ state.

To obtain a quantitative measure of the difference between the irregular and regular evolution, we calculate the overlap of number distributions. The irregular evolution of the number distribution in the chaotic region (see the lower panel of Fig. 9) can be expected to lead to a loss of overlap between number distributions that are initially similar, or between number distributions that are initially identical but which evolve under slightly different Hamiltonians.

We consider the overlap of number distributions, rather than the overlap of the entire state Graß et al. (2013), since (a) number distributions are much more experimentally accessible, and (b) in a theoretical calculation it may be feasible to calculate the number distribution even when the full state is intractable. In particular, we benchmark the use of binned truncated Wigner trajectories Lewis-Swan et al. (2016) to calculate the number distribution approximately, since this method remains tractable for large numbers of particles and modes.

For two states expressed in the number-state basis, $\left|\psi\right>=\sum_{n}c_{n}\left|n\right>$ and $\left|\psi^{\prime}\right>=\sum_{n}d_{n}\left|n\right>$, the state overlap is

[TABLE]

whereas the overlap of number distributions $P_{n}=\left|c_{n}\right|^{2}$ and $\tilde{P}_{n}=\left|d_{n}\right|^{2}$ is given by

[TABLE]

The triangle inequality holds that $B(\psi,\psi^{\prime})\geq O(\psi,\psi^{\prime})$. Although the number distribution overlap is not equivalent to state overlap, it provides an upper bound and can function as a practical alternative.

To obtain a quantum analogue of the maximal Lyaponuv exponent, we calculate the Bhattacharyya statistical distance Bhattacharyya (1943):

[TABLE]

developed as a means of characterising the divergence between two probability distributions and which takes the minimum value of zero when the number distributions are identical.

Following Ref. Lewis-Swan et al. (2016), we construct the approximate number distribution for site 1 by calculating the proportion $P_{n}$ of truncated Wigner trajectories that fall into bins satisfying

[TABLE]

We can expect this binning approach to be accurate when the Wigner function does not change rapidly over the bin size and for a large number of particles such that the corrections from symmetric ordering can be neglected. Fig. 9 gives a comparison of the exact and reconstructed number distributions for the three representative states after 20 driving periods. In all cases, the binned reconstructions recover the overall extent and the coarse features of the distributions. On the other hand, the fine features of the regular 2 and chaotic states are lost in the reconstructions, which raises the question of how sensitive to chaos the Bhattacharyya distance based on these distributions will be.

We calculate the Bhattacharyya distance between distributions evolved under Hamiltonians that differ slightly in the strength of the interaction: $U$ and $U^{\prime}=(1+p)U$, for $p\ll 1$. Fig. 10 shows the results for the binned trajectories calculation (left) and the the direct state-vector calculation (right). Despite the smoothing effect seen in Fig. 9, the Bhattacharyya distance for the chaotic state is distinctly larger than that of the regular states, so long as a sufficiently small perturbation is chosen. However, the Bhattacharyya distances determined via the truncated Wigner and exact state methods quantitatively differ, likely due to the loss of fine features in the number distribution evident in the truncated Wigner reconstruction shown in Fig. 9.

Moreover, the truncated Wigner Bhattacharyya distance for the chaotic state saturates to about $D_{B}\approx 2.4\times 10^{-3}$ for perturbations of $p=10^{-4},10^{-5}$, behaviour which is not evident for the corresponding exact state Bhattacharyya distances. This saturation is consistent with the statistical uncertainty in the binning process. For an estimate of the saturation limit, consider a uniform number distribution $P_{n}=\frac{1}{N+1}$ and a perturbed distribution $\tilde{P}_{n}=P_{n}+\epsilon$, where $\epsilon$ is a Gaussian random variable with mean zero and standard deviation $\sigma$; the Bhattacharyya distance has analytic value $D_{B}=-\ln{\left(1-\frac{(N+1)^{2}}{8}\sigma^{2}\right)}$ in the limit of large $N$. For $D_{B}\approx 2.4\times 10^{-3}$, this estimate gives $\sigma\approx 1.4\times 10^{-4}$, which is approximately consistent with a Poissonian process sorting $N_{t}=10^{5}$ trajectories into $N_{b}=10^{3}$ bins with mean bin occupation $\bar{p}=10^{2}$ and occupation standard deviation $\sigma_{p}=10$ (Poissonian):

[TABLE]

Choosing $10^{4}$ trajectories causes the Bhattacharyya distance to saturate at around $D_{B}\approx 2.4\times 10^{-2}$, also consistent with a Poissonian process for this reduced number of trajectories. Therefore, we can conclude that the saturation of the Bhattacharyya distance is due to sampling error in the binning process.

V Conclusions

In summary, we have explored the quantum and semiclassical dynamics of the modulated Bose-Hubbard dimer around the point at which the semiclassical model undergoes a bifurcation. In the semiclassical model, the extent of the chaos can be controlled by variation of both the modulation frequency and amplitude. Chaotic behaviour tends to be suppressed at high modulation frequency, consistent with the effective static integrable Hamiltonian that can be derived in this regime.

The diverging semiclassical trajectories in the chaotic region correspond, in the quantum dynamics, to phase-space distributions that spread out over time to fill the entire chaotic sea. For a two-mode system, a broad phase-space distribution on the Bloch sphere, corresponding to a diminished average spin vector, correlates directly with fragmentation. We found, nevertheless, that it was states in the regular region undergoing nonlinear shearing that had the greatest loss of condensate fraction (i.e. greater fragmentation).

Despite the complex behaviour of the underlying phase-space distributions, we found that an ensemble of truncated Wigner trajectories reproduced the behaviour of exact quantum expectation values of chaotic quantum states to a high degree of accuracy.

To obtain a measure for quantum systems that plays a similar role that Lyaponuv exponents do in classical chaos, we considered the divergence of number distributions caused by perturbations to the Hamiltonian, as quantified by the Bhattacharyya distance. We benchmarked the use of binned truncated Wigner trajectories as an approximation to the number distribution that is extensible to large-$N$, many site lattices, where exactly calculating the state overlap becomes intractable. We found that for sufficiently small perturbations to the Hamiltonian, the Bhattacharyya distance reliably distinguished between quantum states localised in the regular region from those in the chaotic region. This is in spite of the lack of quantitative agreement in the values of the Bhattacharyya distance calculated with binned Wigner trajectories compared with those calculated exactly. We attribute this discrepancy to the truncated Wigner approach smoothing over some of the fine features of the true number distribution.

Appendix A Effective Hamiltonian for high-frequency modulation

To derive the effective Hamiltonian (5), we introduce the pseudospin operators $\hat{J}_{x,y,z}$ satisfying commutation relations $[\hat{J}_{j},\hat{J}_{k}]=i\epsilon_{jkl}\hat{J}_{l}$, where $j,k,l\in\{x,y,z\}$ and $\epsilon_{jkl}$ is the Levi-Civita symbol,

[TABLE]

The time-independent $H_{0}$ and time-dependent $H_{1}$ parts of the Hamiltonian (1) can then be written as:

[TABLE]

The effective Hamiltonian Kohler et al. (2005); Gong et al. (2009) under a high-frequency modulation of period $T=2\pi/\omega$ is:

[TABLE]

where we have written $\int_{0}^{t}\!\hat{H}_{1}(t^{\prime})\,\mathrm{d}t^{\prime}=A(t)\hat{J}_{x}$, with $A(t)\equiv\frac{2\mu}{\omega}\sin{(\omega t)}$.

To evaluate the self-interaction term in Eq. (17), we can use the Baker-Hasdorff lemma Sakurai and Napolitano (2014):

[TABLE]

together with the integrals:

[TABLE]

where $\mathcal{J}_{0}(z)$ is the zeroth order Bessel function of the first kind.