Enhancement of many-body quantum interference in chaotic bosonic systems

TL;DR

This paper investigates how discrete symmetries enhance many-body quantum interference effects in chaotic bosonic systems, leading to observable deviations from classical approximations like the TWA, with evidence of dynamical localization.

Contribution

It introduces an augmented TWA that incorporates symmetry-induced quantum interference effects, improving the understanding of quantum phenomena in chaotic bosonic systems.

Findings

Enhanced quantum observables due to symmetry effects

Evidence of dynamical localization in Bose-Hubbard model

Deviations from classical TWA predictions indicating quantum interference

Abstract

Although highly successful, the truncated Wigner approximation (TWA) leaves out many-body quantum interference between mean-field Gross-Pitaevskii solutions as well as other quantum effects, and is therefore essentially classical. Turned around, this implies that if a system's quantum properties deviate from TWA, they must be exhibiting some quantum phenomenon, such as localization, diffraction, or tunneling. Here, we consider in detail a particular interference effect arising from discrete symmetries, which can lead to a significant enhancement of quantum observables with respect to the TWA prediction, and derive an augmented version of the TWA in order to incorporate them. Using the Bose-Hubbard model for illustration, we further show strong evidence for the presence of dynamical localization due to remaining differences between the TWA predictions and quantum results.

Equations

$$\hat{H}=\sum _{l,l^{\prime }=1}^L{H}_{l,l^{\prime }}{\hat{b}}_l^{†}{\hat{b}}_{l^{\prime }}+\frac{1}{2}\sum _{l=1}^L{U}_l{\hat{b}}_l^{†}{\hat{b}}_l^{†}{\hat{b}}_l{\hat{b}}_l$$

$${\hat{q}}_j+i{\hat{p}}_j=\sqrt{2\tilde{\hslash }}{\hat{b}}_j,{\hat{q}}_j-i{\hat{p}}_j=\sqrt{2\tilde{\hslash }}{\hat{b}}_j^{†},$$

$$⟨{\mathbf{q}}_f|e^{-i\tau \hat{H}/\tilde{\hslash }}|{\mathbf{q}}_i⟩=\sum _{\gamma }{A}_{\gamma }({\mathbf{q}}_f,{\mathbf{q}}_i)e^{i{R}_{\gamma }({\mathbf{q}}_f,{\mathbf{q}}_i)/\tilde{\hslash }-i{\kappa }_{\gamma }\pi /2}$$

$$i\tilde{\hslash }\frac{\partial }{\partial \tau }{\psi }_l(\tau )=\sum _{l^{\prime }=1}^L{H}_{l,l^{\prime }}{\psi }_{l^{\prime }}(\tau )+U_l(|{\psi }_l(\tau ){|}^2-1){\psi }_l(\tau )$$

$$⟨\hat{f}⟩(\tau )=\int d{\mathbf{q}}_i^{\prime }d{\mathbf{q}}_i^{\prime \prime }d{\mathbf{q}}_f^{\prime }d{\mathbf{q}}_f^{\prime \prime }{\Phi }_0^{\ast }({\mathbf{q}}_i^{\prime }){\Phi }_0({\mathbf{q}}_i^{\prime \prime })⟨{\mathbf{q}}_f^{\prime }|\hat{f}|{\mathbf{q}}_f^{\prime \prime }⟩$$

$$⟨\hat{f}{⟩}_{\mathrm{diag}}(\tau )=\int \frac{d{\mathbf{Q}}_{i}d{\mathbf{P}}_i}{(2\pi \tilde{\hslash }{)}^L}[{\rho }_0{]}_W({\mathbf{Q}}_i,{\mathbf{P}}_i)[f{]}_W({\mathbf{Q}}_f,{\mathbf{P}}_f),$$

$$[O{]}_W(\mathbf{Q},\mathbf{P})\equiv \int d\delta \mathbf{q}{e}^{(i/\tilde{\hslash })\mathbf{P}\cdot \delta \mathbf{q}}⟨\mathbf{Q}+\frac{\delta \mathbf{q}}{2}|\hat{O}|\mathbf{Q}-\frac{\delta \mathbf{q}}{2}⟩.$$

$${\Phi }_0(\mathbf{q})=\frac{1}{(\pi \tilde{\hslash }{)}^{\frac{L}{4}}}\exp \left[-\frac{(\mathbf{q}-{\mathbf{q}}^{(0)}{)}^2}{2\tilde{\hslash }}+\frac{i}{\tilde{\hslash }}{\mathbf{p}}^{(0)}\cdot (\mathbf{q}-{\mathbf{q}}^{(0)})\right].$$

$$A_{\gamma }({\mathbf{q}}_f,{\mathbf{q}}_i)=\frac{e^{i\pi \alpha /4}}{(2\pi \tilde{\hslash }{)}^{L/2}}{\left|\det \left(-\frac{{\partial }^2{R}_{\gamma }}{\partial q_f^l\partial q_i^{l^{\prime }}}({\mathbf{q}}_f,{\mathbf{q}}_i)\right)\right|}^{1/2}$$

$$⟨\hat{f}⟩(t{)}_{\mathrm{diag}}=\int d{\mathbf{Q}}_{i}d{\mathbf{Q}}_{f}d{\delta \mathbf{q}}_{f}d{\delta \mathbf{q}}_i⟨{\mathbf{q}}_i^{\prime \prime }|{\hat{\rho }}_0|{\mathbf{q}}_i^{\prime }⟩⟨{\mathbf{q}}_f^{\prime }|\hat{f}|{\mathbf{q}}_f^{\prime \prime }⟩.$$

Full text

spin echo, many body, semiclassical approximations, fermions, interactions, interference, spin-orbit coupling

For an ultracold bosonic system sufficiently isolated from its environment, a host of fascinating many-body quantum phenomena can be explored. Effects such as various aspects of quantum tunneling Hensinger et al. (2001); Fölling et al. (2007), quantum revivals of coherent states that have initially dispersed Greiner et al. (2002), and the more subtle coherent backscattering in Fock space Engl et al. (2014, 2018), are all significant examples.

In the context of single particle or few-body systems, some phenomena have their most striking manifestation in a system’s time evolution, such as the aforementioned revivals Yeazell et al. (1990); Robinett (2004), but others, such as localization in its various forms, also have unique signatures in the time domain Fishman et al. (1982); Shepelyansky (1983); Heller (1984); Ketzmerick et al. (1992); Bohigas et al. (1993). The onsets of these signatures typically begin beyond an Ehrenfest time scale, $\tau_{E}$, Ehrenfest (1927); Berman and Zaslavsky (1978) in which even the most localized initial states must disperse. For bounded, strongly chaotic dynamical systems, this time scale is logarithmically-short-in-$\hbar$.

The high density limit of ultracold bosonic many-body systems has a deep formal similarity to the semiclassical limit of few-body systems, with the inverse filling factor playing the mathematical role of $\hbar$, and the classical mean-field solutions (of the Gross-Pitaevskii (GP) equation or its discrete version) the role of classical trajectories Engl et al. (2014). A complete semiclassical theory includes interference, diffraction, and tunneling, and thus, generally speaking, quantum phenomena identified in few-body systems are expected to have their analogs in many-body systems. This, in particular, implies a breakdown beyond $\tau_{E}$ of the “naive” classical mean-field approximation. It is far from evident though how such post-Ehrenfest processes can be employed in practice for exploring genuinely quantum effects associated with bosonic many-body systems, such as interferences in Fock space. Indeed, although these interferences are supposed to manifest themselves sensitively in the form of rapid oscillations within the quantum many-body wavefunction’s time evolution, the observable impacts of these rapid oscillations typically are effectively washed out if the detailed information contained in the many-body wavefunctions is projected onto a generic single-particle operator’s expectation value.

This actuality, discussed below in greater detail, lies at the heart of the success of the essentially classical phase-space method known as the truncated Wigner approximation (TWA) Steel et al. (1998); Sinatra et al. (2002); Polkovnikov (2010). In practice, TWA amounts to performing a Monte-Carlo sampling of the time evolution of a quantum many-body state in terms of GP trajectories, i.e., classical fields that evolve according to a GP equation, whose initial values are randomly chosen such that they correctly sample the phase space distribution of the initial quantum state under consideration Steel et al. (1998); Sinatra et al. (2002). The possible occurrence of quantum interference between those GP trajectories is completely neglected within the TWA (as are diffractive or classically prohibited trajectories), more or less as though the systems were weakly, but just sufficiently, coupled to a decohering environment.

As long as some effective time average is performed when comparing with experimental data, there is a general expectation that the above reasoning remains valid also for the expectation values of more sophisticated many-body observables, such as the detection probability of a given Fock state with respect to a single-particle basis. Conversely, a significant deviation of a system’s time-averaged behavior with respect to the TWA prediction is a sensitive indicator of some special many-body quantum phenomenon.

Ahead, coherent state survival probabilities, i.e. the absolute square of an evolved coherent state’s overlap with it’s initial self, are considered as an especially interesting class of measures. The coherent states have minimum uncertainty and are initially the most classical possible Glauber (1963). Furthermore, the survival probability, with some time averaging that eliminates the “generic” rapid interference oscillations contains a great deal of information about the various other surviving and more robust forms of quantum effects. Discrete symmetries, which have a long history, e.g. from Bragg scattering Martin et al. (1988), molecular spectroscopy Harris and Bertolucci (1989), electronic transport Baranger et al. (1991), to topological materials Ando (2013), lead to such surviving constructive interferences. Here we show how to incorporate those effects into the TWA, and how such symmetry-related interferences emerge after $\tau_{E}$. We also argue that the further deviations after time averaging between the TWA and certain Bose-Hubbard model cases shown ahead are signatures of dynamical localization due to quantum interference in a many-body context.

Clearly, for the survival probability an excess of constructive interference for any reason would lead to the TWA underestimating the quantum results, and conversely, an excess of destructive interference would lead to the TWA overestimating them. It is not immediately obvious how to incorporate dynamical effects into the TWA such as scarring Heller (1984), dynamical localization in systems with classical transport barriers Bohigas et al. (1993) or diffusive dynamics Fishman et al. (1982), or tunneling, although such effects would create telltale signals in the behavioral differences between the TWA and quantum results Ullmo et al. (2020). However, it is possible to derive an enhanced version of the TWA accounting for discrete symmetry quantum interference effects.

Indeed, as illustrated in Fig. 1, if both the Hamiltonian and the initial state of the survival probability are symmetric individually with respect to some group of symmetry operations (e.g., reflection, permutation, discrete rotation), then the usual TWA calculation underestimates the mean value of such an observable beyond $\tau_{E}$ by a time-dependent mean symmetry factor $\langle g\rangle$. The unaccounted for enhancement is due to systematic constructive many-body interferences that arise between the contributions of symmetry-related GP trajectories. This enhancement is similar to phenomena found in a single particle physics context Baranger et al. (1991); O’Connor and Tomsovic (1991). However, with increasing numbers of degrees of freedom comes the possibility of much larger discrete symmetry groups where this factor can become arbitrarily large.

The TWA can be obtained from a semiclassical van Vleck-Gutzwiller description of the time evolution of ultracold bosonic system by applying the diagonal approximation Berry (1985); Sun et al. (1998); Dittrich et al. (2006); Dujardin et al. (2015); Tomsovic et al. (2018). To show this, we assume that the time evolution can be well represented within a finite discrete one-body basis containing $L$ single-particle wavefunctions. In the case of a Bose-Einstein condensate that is confined within an optical lattice, those single-particle wavefunctions would most naturally be given by the Wannier orbitals associated with the individual lattice sites. The quantum system is generally described by a Bose-Hubbard Hamiltonian,

[TABLE]

with $H_{l,l^{\prime}}=H_{l^{\prime},l}^{*}$ the one-body Hamiltonian matrix elements and $U_{l}$ the on-site interaction energies.

We introduce the quadrature operators $(\hat{q}_{j},\hat{p}_{j})$ defined as

[TABLE]

where $\tilde{\hbar}$ equals the mean filling factor inverse, i.e. ratio of the site number to mean total particle number, $L/N$. Following Ref. Dujardin et al. (2015), we represent the semiclassical propagator in the $|\mathbf{q}\rangle$ basis associated with $\mathbf{\hat{q}}\equiv(\hat{q}_{1},\ldots,\hat{q}_{L})$. For a large total particle number, the time evolution operator matrix elements can be expressed within the quadrature basis as a sum over mean-field trajectories Littlejohn (1992),

[TABLE]

which go from $\mathbf{q}_{i}$ to $\mathbf{q}_{f}$ in the scaled time $\tau=\tilde{\hbar}t/\hbar$. For the Hamiltonian of Eq. (1) such a trajectory would correspond to a solution of the GP equation

[TABLE]

satisfying the boundary conditions $\mathrm{Re}[\psi^{l}(0)]=q^{l}_{i}/\sqrt{2\tilde{\hbar}}$ and $\mathrm{Re}[\psi^{l}(\tau)]=q^{l}_{f}/\sqrt{2\tilde{\hbar}}$. $R_{\gamma}(\mathbf{q}_{f},\mathbf{q}_{i})$ represents its associated principal function or action integral. The integer Morse index $\kappa_{\gamma}$ Littlejohn (1992) counts the number of conjugate points along the trajectory 111See supplemental material for further details..

Any generic single- or many-particle observable in this bosonic system can be expressed in terms of the quadrature operators as $\hat{f}\equiv f(\hat{\mathbf{q}},\hat{\mathbf{p}})$ yielding, e.g., $\hat{b}_{l}^{\dagger}\hat{b}_{l}=(\hat{q}_{l}^{2}+\hat{p}_{l}^{2}-\tilde{\hbar})/2\tilde{\hbar}$ for the occupancy on site $l$. Given an initial many-body state $|\Phi_{0}\rangle$, the time evolution $\langle\hat{f}\rangle(\tau)\equiv{\rm Tr}[\hat{f}\hat{\rho}_{\tau}]$ of the mean value of such an observable, with $\hat{\rho}_{\tau}\equiv|\Phi_{\tau}\rangle\langle\Phi_{\tau}|$ and $|\Phi_{\tau}\rangle\equiv e^{-i\tau\hat{H}/\tilde{\hbar}}|\Phi_{0}\rangle$, is then expressed as

[TABLE]

where $\Phi_{0}(\mathbf{q})$ represents the initial quantum many-body wavefunction in the $q$-quadrature basis.

For an averaging process that suppresses oscillating terms, the main contributions to this integral are expected to be given by the diagonal approximation. This implies two crucial assumptions: (i) that in the double sum over trajectories, all terms for which the two action integrals $R_{\gamma^{\prime}}$ and $R_{\gamma^{\prime\prime}}$ correspond to different orbits cancel out; and (ii) only short chords Ozorio de Almeida (1998), i.e. points such that $\mathbf{q}_{i}^{\prime}\simeq\mathbf{q}_{i}^{\prime\prime}$ and $\mathbf{q}^{\prime}_{f}\simeq\mathbf{q}_{f}^{\prime\prime}$ are going to contribute significantly to the integrals. If so, one may transform the double sum in Eq. (Enhancement of many-body quantum interference in chaotic bosonic systems: the role of symmetry and dynamics) into a single sum over mean field trajectories and expand in the small parameter associated with the chord. This leads to

[TABLE]

where $(\mathbf{Q}_{f},\mathbf{P}_{f})$ has to be understood as the final point in phase space at time $\tau$ of the trajectory initiated at $(\mathbf{Q}_{i},\mathbf{P}_{i}$), and with the Wigner transform of an operator $\hat{O}$ defined as

[TABLE]

Provided $f$ is a rather smooth and well-behaved function of $\mathbf{q}$ and $\mathbf{p}$, $[f]_{W}(\mathbf{Q}_{f},\mathbf{P}_{f})$ can, to within small $\tilde{\hbar}$ corrections, be approximated by $f(\mathbf{Q}_{f},\mathbf{P}_{f})$, and the integral in Eq. (6) can be evaluated by a Monte-Carlo method. This becomes the TWA.

The above reasoning would be valid if either the expectation value itself, or some further temporal or configuration averages, removes the contributions of Eq. (Enhancement of many-body quantum interference in chaotic bosonic systems: the role of symmetry and dynamics)’s off-diagonal terms. However, if there exists a discrete symmetry, more care must be exercised. Consider an initial condensate wavefunction that is symmetric with respect to a given parity exhibited by the Hamiltonian and a many-body observable $\hat{f}$. Within the phase space there is a subspace of points that are their own parity partners. Denote it the symmetry subspace. At long times nearly every trajectory $\gamma$ that significantly contributes to the semiclassical expression, Eq. (Enhancement of many-body quantum interference in chaotic bosonic systems: the role of symmetry and dynamics), for the mean value of this observable leaves the neighborhood of the symmetry subspace and has a symmetry-related partner trajectory $\tilde{\gamma}$, which is obtained by applying the parity operator onto $\gamma$ and which exhibits the same action integral as $\gamma$. In particular, partner trajectories satisfy $R_{\tilde{\gamma}}(\tilde{\mathbf{q}}_{{f}},\tilde{\mathbf{q}}_{{i}})=R_{\gamma}(\mathbf{q}_{{f}},\mathbf{q}_{{i}})$ and thus also $A_{\tilde{\gamma}}(\tilde{\mathbf{q}}_{{f}},\tilde{\mathbf{q}}_{{i}})=A_{\gamma}(\mathbf{q}_{{f}},\mathbf{q}_{{i}})$. The contributions of those two trajectories therefore interfere constructively within Eq. (Enhancement of many-body quantum interference in chaotic bosonic systems: the role of symmetry and dynamics). As a consequence, Eq. (6), which entirely neglects those interferences, underestimates the true expectation value of $\hat{f}$ by a factor two.

To account for this effect in a quantitatively correct manner within TWA, each trajectory must be associated with its particular symmetry factor, $g_{\gamma}$, that correctly counts the number of other symmetry-related trajectories with which constructive interference will arise. A naive way to approach that problem would be to perform the Monte-Carlo calculation implied by Eq. (6) and multiply each contribution by its number of existing distinct symmetry related orbits. However, since the number of orbits that are symmetric under an element of the symmetry group is of measure zero, this would lead to multiplying the TWA result by a global factor $g_{\rm max}$ (the number of symmetry group elements), which entirely misses the transient shift from unity to maximum at relatively short dynamical times, and misses the fact that due to coherent state spread, the long time $\langle g\rangle$ can saturate below $g_{\rm max}$; see Fig. 1a where $g_{\rm max}=4$.

A less intuitive and more precise picture for why this approach cannot succeed is that the correct enhancement factor $g_{\gamma}$ is not the symmetry of the orbit itself, but that of its family of orbits defined by the neighborhood of a “saddle trajectory” Tomsovic (2018). This aspect appears when some of the integrals of Eq. (Enhancement of many-body quantum interference in chaotic bosonic systems: the role of symmetry and dynamics) are performed within the stationary phase approximation (hence, “saddle trajectory”), and the contributions of orbit families (saddle trajectory neighborhoods) are naturally grouped together. This grouping ensures that orbits within a family are similar for their entire evolution, and in particular, are well represented approximately through the use of the saddle trajectory’s stability matrix. If an orbit and a symmetric version are within the same family, there is no enhancement. If they are associated with well separated saddles, which cannot occur in the dynamics for times shorter than $\tau_{E}$, then there is an enhancement. Note that the stationary phase approximation involves locating complex saddle trajectories Huber and Heller (1987); Huber et al. (1988); Pal et al. (2016); Baranger et al. (2001); Tomsovic et al. (2018); Tomsovic (2018) in order to evaluate the mean value of observables, and lacks the simplicity of the “initial value” formulation of the TWA, Eq. (6).

Nevertheless, the conventional implementation of the TWA can be amended to incorporate the impact of discrete symmetries in an automated manner. First, note that by the saddle logic above, the “symmetric” trajectories with respect to a given symmetry operation that do not lead to enhancement must remain during their entire evolution in the neighborhood of the relevant symmetry subspace. Those non-symmetric ones leading to an enhancement factor do not. Hence, in order to discriminate between those possibilities, it is sufficient to introduce an $\tilde{\hbar}$ distance scale $d$ in phase space through which the notion of “closeness” to a symmetry subspace can be properly defined. A trajectory that is computed within the framework of Eq. (6) is symmetric with respect to a given symmetry subspace if its distance to that subspace remains below $d$ for its entire evolution time. Otherwise, it is non-symmetric with respect to that subspace, and in that case its contribution is multiplied by the associated enhancement factor. Up to $\tau_{E}$, trajectories cannot leave a symmetric subspace and return. The mean value of $f$ is then provided by symmetric trajectories, i.e., pairs, both of whose members are included within the short chord approximation: their $\gamma$ does not give rise to a distinct $\tilde{\gamma}$ ($\gamma=\tilde{\gamma}$). Thus, the the dynamics need time to explore phase space sufficiently to manifest the symmetry’s existence, and this leads to time scales associated with increasing multiplicities, depending on the symmetry group.

The “symmetry-enhanced” curves in Fig. 1 have been constructed following this approach. It shows a homogeneous $1D$ Bose-Hubbard chain with $L=4$ sites and periodic boundary conditions. This system is modeled by Eq. (1) with the specific choices $U_{l}=U>0$ for all $l$ and $H_{ll^{\prime}}=-J<0$ if $l^{\prime}=l\pm 1$ (mod $L$) and zero otherwise. Furthermore, consider an initial coherent state for a perfect BEC centered about the field amplitudes $\psi_{l}^{(0)}=(q_{l}^{(0)}+ip_{l}^{(0)})/\sqrt{2\tilde{\hbar}}$ 222The breaking of $U(1)$ symmetry that is implied by this ansatz does not play a significant role for most commonly considered observables, provided the BEC contains a reasonably large number $N\gg 1$ of particles Lieb et al. (2007), which is equivalent to a wave packet up to a global phase Tomsovic et al. (2018)

[TABLE]

More specifically, consider the coherent state density waves: (a) $(\sqrt{9},1,\sqrt{9},1)$ and (b) $(\sqrt{18},\sqrt{2},\sqrt{9},\sqrt{2})$ as initial states, and the many-body observable of interest given by $\hat{f}=|\Phi_{0}\rangle\langle\Phi_{0}|=\hat{\rho}_{0}$, whose expectation value is the survival probability, $|\langle\Phi_{0}|\Phi_{\tau}\rangle|^{2}$, after a given evolution time $\tau$, where we perform a time average rem in order to filter out rapid oscillations arising due to the generic quantum interferences. The symmetry group in this example is larger than the one associated with simple parity. As a consequence, the enhancement factor, which depends on the group, is larger than two. Indeed, for a Bose-Hubbard ring with equal on-site interaction energies (all $U_{l}=U$) and only equal nearest neighbor one-body terms, the full discrete group includes cyclic permutations and a clockwise-counterclockwise equivalence. The maximum possible enhancement factor for such a system, before accounting for the symmetry of the initial state and observable, is twice the number of sites, $2L$. However, the density wave partially breaks the full system dynamical symmetry as only even cyclic permutations leave it invariant. Thus, the maximum symmetry enhancement factor is reduced by a factor two to $L$, which equals $4$ in our figures.

The curve thus obtained accounting for $g_{\gamma}$ (noisy red line) is seen in Fig. 1 to perfectly follow the smoothed exact quantum result. Evidently, this approach can straightforwardly be generalized to account for the presence of larger symmetry groups in the Hamiltonian and the initial state. In calculations not shown with $L=6$, the enhanced TWA saturated at a factor $\langle g\rangle=6$, and the agreement with the smoothed exact quantum result is excellent. Note that this factor could be much larger. For the fully connected Bose-Hubbard model studied in Ref. Pizzi et al. (2019), the maximum enhancement factor for a density wave survival probability would be $(L/2)!^{2}$.

In a second example shown in Fig. 2, a different initial coherent state is chosen, other parameters the same as before, for which just beyond $\tau_{E}$ there is rough agreement between the symmetry enhanced TWA and the smoothed quantum survival probability. However, with increasing time there is a decrease in the TWA indicating that new regions of phase space are opening up with time, and yet the quantum average does not follow this decrease. This demonstrates that the system is undergoing additional quantum interference effects, which in this case strongly suggests a dynamical localization process. One for which there are time scales beyond the Ehrenfest time in which the classical trajectories are continuing to access new regions of phase space that the quantum system cannot.

In summary, the very successful TWA approximation can be extended to account quantitatively for basic constructive interference effects enforced by the existence of discrete symmetries, which can involve large enhancement factors. In the time domain, discrete symmetries reveal themselves over time and depending on the dynamics incorporate multiple time scales. For example, the $L=8$ ring has subspaces with factors, $1,2,4,8$. For some initial condensate coherent states and values of $NU/LJ$, the Lyapunov exponent in the direction of the double degeneracy symmetry subspace is much greater than that for $4$ or $8$. The factor $2$ then enters the dynamics well before the higher degeneracy factors.

In addition, deviations between the TWA and quantum behaviors of the survival probability indicate the presence of some genuine quantum effect such as localization, diffraction, or tunneling. In the case shown here, the initial agreement followed by a divergence appears to be due to a quantum interference effect leading to dynamical localization. This is distinctly different from a classical localization in which the GP solutions are trapped or fail to explore phase space regions common in mixed phase space systems, for example, the macroscopic self-trapping discussed in Ref. Pizzi et al. (2019). More work is needed to separate signals of the various forms of localization from each other or from tunneling and diffraction effects, and is left for the future.

Supplementary material

In this supplementary material we provide some of the details leading to Eq. (6). First, note that the prefactor can be written as

[TABLE]

with $\alpha$ the index of inertia of $-\partial\mathbf{q}_{i}/\partial\mathbf{p}_{i}$ Littlejohn (1992). Next, starting from Eq. (5), assume that the contributions of non-identical pairs of orbits cancel out because of the time average, and that only short chords, i.e. points such that $\mathbf{q}_{i}^{\prime}\simeq\mathbf{q}_{i}^{\prime\prime}$ and $\mathbf{q}^{\prime}_{f}\simeq\mathbf{q}_{f}^{\prime\prime}$ contribute significantly. In such a case, the variables $\mathbf{Q}_{f,i}\equiv(\mathbf{q}_{f,i}^{\prime\prime}+\mathbf{q}_{f,i}^{\prime})/2$ and $\mathbf{\delta q}_{f,i}\equiv(\mathbf{q}_{f,i}^{\prime\prime}-\mathbf{q}_{f,i}^{\prime})$ can be introduced. The key is to expand in the “small” variables $\mathbf{\delta q}_{f,i}$, which means more specifically keeping only the zero’th order term for the smooth function $A_{\gamma}$ but expanding the action $R_{\gamma}$ to first order. With $\gamma^{\prime}=\gamma^{\prime\prime}=\gamma$, and using the property $\mathbf{p}^{(\gamma)}_{f}=\partial R_{\gamma}/\partial\mathbf{q}_{f}$; $\mathbf{p}^{(\gamma)}_{i}=-\partial R_{\gamma}/\partial\mathbf{q}_{i}$ (with $\mathbf{p}^{(\gamma)}_{i}$ and $\mathbf{p}^{(\gamma)}_{f}$ the initial and final “momenta” of the trajectory $\gamma$), gives

[TABLE]

with $\hat{\rho}_{0}\equiv|\Phi_{0}\rangle\langle\Phi_{0}|$ the initial density. Note that $-{\partial^{2}R_{\gamma}}/{\partial\mathbf{Q}_{f}\partial\mathbf{Q}_{i}}={\partial\mathbf{P}_{i}}/{\partial\mathbf{Q}_{f}}_{\gamma}$. Therefore, the determinant in Eq. (9) is just the Jacobian of the transformation from the final “position” $\mathbf{Q}_{f}$ to the initial “momentum” $\mathbf{P}_{i}$. Thus, $\sum_{\gamma}\int d\mathbf{Q}_{f}|A_{\gamma}|^{2}\mapsto(2\pi\tilde{\hbar})^{L}\int d\mathbf{P}_{i}$ leading to Eq. (6).