Phase-space mixing in dynamically unstable, integrable few-mode quantum systems

TL;DR

This paper investigates phase-space mixing and relaxation dynamics in quantum few-mode systems that are classically integrable but become dynamically unstable after a quench, using analytical methods and numerical simulations.

Contribution

It provides analytical expressions for observable dynamics and long-time limits in unstable, integrable quantum systems, validated by numerical simulations.

Findings

Long-time expectation value deviation scales as 1/ln(N).

Observable relaxation is a Gaussian-damped oscillation.

Results confirmed by numerical TWA simulations.

Abstract

Quenches in isolated quantum systems are currently a subject of intense study. Here, we consider quantum few-mode systems that are integrable in their classical mean-field limit and become dynamically unstable after a quench of a system parameter. Specifically, we study a Bose-Einstein condensate (BEC) in a double-well potential and an antiferromagnetic spinor BEC constrained to a single spatial mode. We study the time dynamics after the quench within the truncated Wigner approximation (TWA) and find that system relaxes to a steady state due to phase-space mixing. Using the action-angle formalism and a pendulum as an illustration, we derive general analytical expressions for the time evolution of expectation values of observables and their long-time limits. We find that the deviation of the long-time expectation value from its classical value scales as $1/O(\ln N )$, where $N$ is the number of atoms in the condensate. Furthermore, the relaxation of an observable to its steady state value is a damped oscillation and the damping is Gaussian in time. We confirm our results with numerical TWA simulations.

Equations

$$⟨\mathcal{O}(t)⟩$$

$$={\int }_{\Omega }d{\mathbf{x}}_{0}d{\mathbf{p}}_0\mathcal{O}(t)F_0({\mathbf{x}}_0,{\mathbf{p}}_0),$$

$${\dot{I}}_i=-\frac{\partial H(\mathbf{I})}{\partial {\phi }_i}=0,{\dot{\phi }}_i=\frac{\partial H(\mathbf{I})}{\partial I_i}\equiv {\omega }_i(\mathbf{I}),$$

$$\mathbf{\phi }(t)=\mathbf{\omega }(\mathbf{I})t+{\mathbf{\phi }}_0,$$

$$⟨\mathcal{O}(t)⟩=\sum _R{\int }_{{\mathcal{I}}_R}d\mathbf{I}{\int }_{\mathcal{J}}\frac{d\mathbf{\phi }}{(2\pi {)}^n}{f}_R(\mathbf{I},\mathbf{\phi },t){\mathcal{O}}_R(\mathbf{I},\mathbf{\phi })$$

$$=\sum _R{\int }_{{\mathcal{I}}_R}d\mathbf{I}{\int }_{\mathcal{J}}\frac{d{\mathbf{\phi }}_0}{(2\pi {)}^n}{f}_{0,R}(\mathbf{I},{\mathbf{\phi }}_0){\mathcal{O}}_R(\mathbf{I},\mathbf{\phi }(t)),$$

$${\bar{f}}_R(\mathbf{I})\equiv {\int }_{\mathcal{J}}\frac{d\mathbf{\phi }}{(2\pi {)}^n}{f}_R(\mathbf{I},\mathbf{\phi },t)={\int }_{\mathcal{J}}\frac{d\mathbf{\phi }}{(2\pi {)}^n}{f}_{0,R}(\mathbf{I},\mathbf{\phi }).$$

$$\underset{t\to \infty }{\lim }⟨\mathcal{O}(t)⟩=\sum _R{\int }_{\mathcal{I}}d\mathbf{I}{\bar{f}}_R(\mathbf{I}){\int }_{\mathcal{J}}\frac{d\mathbf{\phi }}{(2\pi {)}^n}{\mathcal{O}}_R(\mathbf{I},\mathbf{\phi }).$$

$$H_{\mathrm{pend}}=\frac{p^2}{2}+1-\cos \theta ,$$

$${\theta }_{\mathrm{S}\pm }(t)=\pm 2\arcsin (\tanh t),p_{\mathrm{S}\pm }(t)=\pm 2\mathrm{sech}(t).$$

$$p_R(t)\sim \sum _{n=-\infty }^{\infty }\left[\sum _{s=\{S\pm \}}{\chi }_R(s)p_s\left(t-t_{0,R}(s)-n{T}_R\right)\right],$$

$${\mathcal{O}}_R(t)\sim {\mathcal{O}}_{\mathrm{sp}}+\sum _{n=-\infty }^{\infty }\sum _s{\chi }_R(s){\mathcal{D}}_s(t-t_{0,R}(s)-n{T}_R).$$

$${\int }_0^{2\pi }\frac{d{\phi }_1}{2\pi }{\mathcal{O}}_R(\mathbf{I},{\phi }_1)$$

$$\times \sum _{n=-\infty }^{\infty }{\int }_0^{T_R}dt{\mathcal{D}}_s(t-t_{0,R}(s)-n{T}_R).$$

$${\int }_0^{2\pi }\frac{d{\phi }_1}{2\pi }{\mathcal{O}}_R(\mathbf{I},{\phi }_1)$$

$$\times {\int }_{-\infty }^{\infty }dt{\mathcal{D}}_s(t).$$

$$\underset{t\to \infty }{\lim }⟨\mathcal{O}(t)⟩\sim {\mathcal{O}}_{\mathrm{sp}}+\sum _R\frac{{⟨{\omega }_1⟩}_R}{2\pi }\left[\sum _s{\chi }_R(s){\int }_{-\infty }^{\infty }dt{\mathcal{D}}_s(t)\right],$$

$$F_0(\theta ,p)=\frac{1}{2\pi d^2}{e}^{-(x^2+p^2)/(2d^2)},$$

$$\underset{t\to \infty }{\lim }⟨p^2(t)⟩$$

$$⟨\varpi ⟩\equiv {\int }_0^{\infty }d\varpi \varpi \mathcal{F}(\varpi ),$$

$$\mathcal{F}(z)$$

$$\mu \equiv ⟨\varpi ⟩$$

$$\sigma$$

$${\mathcal{O}}_R(\mathbf{I},{\phi }_1)=\sum _{m=-\infty }^{\infty }{\Theta }_R(\mathbf{I};m)e^{im{\phi }_1},$$

$${\Theta }_R(\mathbf{I};m)={\int }_0^{2\pi }\frac{d{\phi }_1}{2\pi }{\mathcal{O}}_R(\mathbf{I},{\phi }_1)e^{-im{\phi }_1}.$$

$${\Theta }_R(\mathbf{I};m)$$

$$\sum _{n=-\infty }^{\infty }{\int }_{-t_{0,R}(s)}^{T_R-t_{0,R}(s)}\frac{d\tau }{2\pi }{\omega }_1{\mathcal{D}}_s(\tau -n{T}_R)e^{-im{\omega }_1\tau },$$

$${\Theta }_R(\mathbf{I};m)\sim {\mathcal{O}}_{\mathrm{sp}}{\delta }_{m0}$$

$$⟨\mathcal{O}(t)⟩$$

$$+\sum _{m=-\infty }^{\infty }(-1{)}^m{\int }_0^{\infty }d\varpi \mathcal{F}(\varpi )\varpi {\mathfrak{D}}_{S+}(m\varpi )e^{im\varpi t},$$

$$⟨\mathcal{O}(t)⟩\sim {\mathcal{O}}_{\mathrm{sp}}+\sum _{m=-\infty }^{\infty }(-1{)}^m\mu {\mathfrak{D}}_{S+}(m\mu )e^{im\mu t-m^2{\sigma }^2{t}^2/2}.$$

Full text

I Introduction

The advent of precise experimental control in ultracold atomic systems has motivated theoretical study in non-equilibrium dynamics in isolated quantum systems Bloch et al. (2008). For generic Hamiltonian systems the expectation value of a local observable at long times after a quench, a sudden change in a control parameter, is described by a Gibbs ensemble Dziarmaga (2010); D’Alessio et al. (2016). However, for integrable systems, a special class of Hamiltonian systems, the long-time behavior is instead believed to be described by a generalized Gibbs ensemble D’Alessio et al. (2016). This important role of integrability on the time dynamics has been demonstrated experimentally Kinoshita et al. (2006); Langen et al. (2015). Integrable systems are of much theoretical interest as they are amenable to exact analytic treatment. A classical integrable system can be solved using action-angle variables Arnold (1997), while a quantum integrable system is solvable by the Bethe ansatz Korepin et al. (1997).

A mean-field approximation can be applied to a bosonic system with a macroscopically-occupied mode. The time dynamics of the system is then governed by a classical Hamiltonian and described by classical trajectories in its phase space. For a weakly interacting Bose-Einstein condensate (BEC) this classical trajectory is a solution of the time-dependent Gross-Pitaevskii equation for the order parameter with continuous spatial degrees of freedom Pethick and Smith (2008). In certain cases, it is sufficient to describe a bosonic system with just a few degrees of freedom. Some examples are a BEC in a double-well potential Smerzi et al. (1997), a spin-1 spinor BEC within the single-mode approximation (SMA) Law et al. (1998); Zhang et al. (2005) and a few-site Bose-Hubbard model with a large occupation per site Mossmann and Jung (2006); Itin and Schmelcher (2011); Satija et al. (2013).

A bosonic system becomes dynamically unstable when it is prepared by a quench at a saddle point in its phase space. Dynamical instabilities have been predicted for vortices in trapped BECs Pu et al. (1999); Skryabin (2000); Kawaguchi and Ohmi (2004), superfluid flow of BECs in optical lattices Wu and Niu (2001); Menotti et al. (2003) and BECs in cavities Moore et al. (1999). These predictions have been experimentally observed Shin et al. (2004); Fallani et al. (2004); Cristiani et al. (2004); Ferris et al. (2008); Schmidt et al. (2014). The instability is also used as an experimental route for the generation of squeezed states Estève et al. (2008); Leslie et al. (2009); Klempt et al. (2010); Hamley et al. (2012). A mean-field description is then insufficient and quantum fluctuations need to be included. Quantum corrections can be (partially) included by using the truncated Wigner approximation (TWA) Sinatra et al. (2002); Blakie et al. (2008); Polkovnikov (2010), which models the dynamics of the Wigner distribution in the phase space. The TWA has been used to numerically study the effects of thermal fluctuations on a BEC Sinatra et al. (2002), quenches in spinor condensates Sau et al. (2010); Barnett et al. (2011) and superfluid flow Mathey et al. (2014).

In this paper we analytically study the time dynamics of two integrable few-mode quantum systems within the truncated Wigner approximation after a quench of a parameter that makes the systems dynamically unstable. Our paper is set up as follows. We introduce dynamical instability in bosonic systems in Sec. II and TWA in Sec. III. We define the integrability of classical Hamiltonians, which govern the mean-field limit of these systems, and introduce action-angle coordinates in Sec. IV. Section V introduces the concept of mixing in phase-space due to time evolution and describes how this mixing leads to relaxation of an observable to a steady-state value. Using the pendulum as an illustrative example, we derive general results for long-time expectation value of an observable in Secs. VI and VII and the time dynamics of relaxation of this expectation value in Sec. VIII. We apply these results to the case of a condensate in a double-well potential (the double-well system) in Sec. IX and a spin-1 BEC described by a single spatial mode in Sec. X. Finally, we conclude in Sec. XI.

II Dynamical instability

The mean-field equations of motion of an isolated quantum bosonic system are equivalent to Hamilton’s equations of motion of a classical system. The mean-field ground state is a stable equilibrium phase-space point, where the classical Hamiltonian has a minimum. On the other hand, a dynamically unstable state corresponds to a saddle point of this Hamiltonian. Such an unstable state can be prepared by starting from a minimum of the initial Hamiltonian and then quenching a system parameter to change this point to a saddle point of the final Hamiltonian. As an example, consider the quantum oscillator $H_{0}=(\hat{p}^{2}+\hat{x}^{2})/2$, where $\hat{x}$ and $\hat{p}$ are the canonical position and momentum operators, respectively. Here, we have set $\hbar$ and the natural frequency of the oscillator to one. Its mean-field ground state is the phase-space point $(x_{c},p_{c})=(0,0)$, where $x_{c}=\left\langle\hat{x}\right\rangle$, $p_{c}=\left\langle\hat{p}\right\rangle$, and $\left\langle\dots\right\rangle$ is the average over a quantum state. We make the state dynamically unstable by suddenly changing to the Hamiltonian $H_{1}=(\hat{p}^{2}-\hat{x}^{2})/2$. Under the mean-field equations of motion a dynamically unstable point is stationary. Thus, $x_{c}(t)=0$ and $p_{c}(t)=0$ holds for all times. In contrast, quantum evolution under $H_{1}$ leads to exponential growth in the unstable mode Pethick and Smith (2008). In fact, following the language of quantum optics, $H_{1}\propto\hat{a}\hat{a}+\hat{a}^{\dagger}\hat{a}^{\dagger}$ leads to single-mode squeezing, where $\hat{a}(\hat{a}^{\dagger})=(\hat{x}\pm i\hat{p})/\sqrt{2}$ is the annihilation (creation) operator of the mode.

III Truncated Wigner Approximation

The time evolution of a dynamically-unstable system can be studied using the truncated Wigner approximation (TWA) Sinatra et al. (2002). It incorporates the leading order quantum corrections to the mean-field equations of motion Polkovnikov (2003). In the TWA a Wigner distribution function $F(\mathbf{x},\mathbf{p},t)$ time evolves under classical Hamilton’s equations, in contrast to the mean-field approximation where the evolution of a single phase-space point $(\mathbf{x}(t),\mathbf{p}(t))$ is studied. Here, $\mathbf{x}=(x_{1},\dots,x_{n})$ and $\mathbf{p}=(p_{1},\dots,p_{n})$ are canonical position and momentum coordinates for a classical mean-field Hamiltonian system with $n$ degrees of freedom. The initial distribution, $F_{0}(\mathbf{x},\mathbf{p})$, is the Wigner transform Case (2008) of the prequench quantum ground state or any approximation thereof.

For an observable $\mathcal{O}(\mathbf{x},\mathbf{p})$, we define its evolution $\mathcal{O}(t)\equiv\mathcal{O}(\mathbf{x}(t),\mathbf{p}(t))$ along a trajectory $(\mathbf{x}(t),\mathbf{p}(t))$ with initial conditions $(\mathbf{x}_{0},\mathbf{p}_{0})$. The expectation value of $\mathcal{O}(t)$ over all trajectories is

[TABLE]

with measures $d\mathbf{x}=dx_{1}\cdots dx_{n}$, $d\mathbf{p}=dp_{1}\cdots dp_{n}$ and the integral is over all phase space $\Omega$. The distribution satisfies $\int_{\Omega}d\mathbf{x}d\mathbf{p}\,F(\mathbf{x},\mathbf{p},t)=1$ for all $t$ in accordance with Liouville’s theorem Arnold (1997).

IV Classical integrable systems

In classical mechanics a Hamiltonian system with $n$ degrees of freedom is integrable if there exist $n$ mutually commuting (with respect to the Poisson bracket) conserved quantities Arnold (1997). Then a trajectory in the $2n$ dimensional phase-space lies on an $n$-dimensional torus. For an integrable system, the coordinates $(\mathbf{x},\mathbf{p})$ can be transformed to canonical coordinates called actions $\mathbf{I}=(I_{1},\dots,I_{n})$ and angles $\bm{\varphi}=(\varphi_{1},\dots,\varphi_{n})$, such that Hamiltonian $H$ is independent of $\bm{\varphi}$. Crucially, $(\mathbf{I},\bm{\varphi})$ and $(\mathbf{I},\bm{\varphi}+2\pi\mathbf{m})$ correspond to the same phase-space point, where $\mathbf{m}=(m_{1},\dots,m_{n})$ is a vector of integers. In these coordinates, the Hamilton’s equations are

[TABLE]

for all $i\in\{1,\dots,n\}$. The frequencies $\omega_{i}(\mathbf{I})$ only depend on $\mathbf{I}$. Hence, the actions are conserved quantities and the time evolution of the angles has the simple form

[TABLE]

where $\bm{\omega}(\mathbf{I})=(\omega_{1}(\mathbf{I}),\dots,\omega_{n}(\mathbf{I}))$ and $\bm{\varphi}(0)=\bm{\varphi}_{0}$.

For our Hamiltonian systems action-angle coordinates are not globally defined. Instead, they are defined on disjoint regions of $\Omega$ by maps from each such region $R$ to $\mathcal{I}_{R}\otimes\mathcal{J}$, where $\mathcal{I}_{R}\subset\mathbb{R}^{n}$ and $\mathcal{J}=[0,2\pi]^{\otimes n}$ are the spaces spanned by the actions and angles, respectively. We then construct distribution functions $f_{R}(\mathbf{I},\bm{\varphi},t)=(2\pi)^{n}F\left(\mathbf{x}(\mathbf{I},\bm{\varphi}),\mathbf{p}(\mathbf{I},\bm{\varphi}),t\right)$ for $(\mathbf{x},\mathbf{p})\in R$ with normalization $\sum_{R}\int_{\mathcal{I}_{R}}d\mathbf{I}\int_{\mathcal{J}}d\bm{\varphi}/(2\pi)^{n}\,f_{R}(\mathbf{I},\bm{\varphi},t)=1$. The latter follows from the fact that canonical transformations have a unit Jacobian. The distribution $f_{R}(\mathbf{I},\bm{\varphi},t)$ is periodic in $\bm{\varphi}$ and evolves as $f_{R}(\mathbf{I},\bm{\varphi},t)=f_{0,R}(\mathbf{I},\bm{\varphi}-\bm{\omega}t)$, where $f_{0,R}(\mathbf{I},\bm{\varphi})=f_{R}(\mathbf{I},\bm{\varphi},0)$ is the initial distribution. Moreover, Eq. 1 becomes

[TABLE]

where $\mathcal{O}_{R}(\mathbf{I},\bm{\varphi})$ is the functional form of the observable in region $R$.

V Phase-space mixing

A distribution function that is initially localized around a phase-space point typically stretches, tangles and disperses over the accessible phase space. This mixing in phase space has been studied in plasma physics Hammett et al. (1992) and astrophysics Tremaine (1999). We illustrate this concept using an anharmonic oscillator. Its Hamiltonian $H=r^{2}/2+\varepsilon r^{4}$ is integrable, where $r^{2}=p^{2}+x^{2}$ and we have set the mass and the natural frequency of the oscillator to unity. In this case the action-angle coordinates are globally defined. The action $I$ is a function of $r$ and the angle $\varphi$ is the polar angle in the $(x,p)$ plane. Points with different $r$ rotate around the origin at different frequencies $\omega(I)$ and the distribution function stretches as shown Fig 1. Eventually, the distribution spreads uniformly and mixes in the compact coordinate $\varphi$, while remaining localized in $r$ and $I$.

For a general integrable system, the frequencies $\bm{\omega}(\mathbf{I})$ depend nontrivially on $\mathbf{I}$. Hence, the distribution will eventually mix in $\bm{\varphi}$. It is important to realize that as the distribution function mixes in phase space fine-scale structures must develop in order to conserve the phase-space volume as required by Liouville’s theorem. For the anharmonic oscillator evolution leads to tightly wound spirals as shown in the third panel of Fig. 1.

Phase-space mixing simplifies the evaluation of the long-time expectation value of an observable. Experimentally-accessible observables are typically smooth functions of the phase-space coordinates. Then the distribution function with its fine-scale structures can be coarsened, i.e., in Eq. 4 we can replace $f_{R}(\mathbf{I},\bm{\varphi},t)$ by the time-independent distribution (Lifshitz and Pitaevskii, 1981, §1)

[TABLE]

Consequently, the expectation value at long times becomes

[TABLE]

Thus, the long-time expectation value of an observable is given by the average over the accessible phase space weighted by $\bar{f}_{R}(\mathbf{I})$.

VI Dynamics near a separatrix

The description of the time evolution of the initially-localized Wigner distribution following dynamical instability for our double-well and spin-1 boson systems with a four- and six-dimensional phase space, respectively, must include a study of separatrices. As we will show in Sec. IX and X their dynamics is controlled by a two-dimensional subspace $\Omega_{\rm 2D}$ spanned by canonical coordinates $x_{1}$ and $p_{1}$. This subspace contains a single saddle point that is connected to itself by one or more trajectories, known as separatrices. In fact, there are two separatrices and one separatrix for the double-well and spin-1 Bose system, respectively. The frequency $\omega_{1}(\mathbf{I})$ associated with a trajectory in $\Omega_{\rm 2D}$ goes to zero as its starting point approaches the saddle point. In fact, near the saddle point $\omega_{1}$ varies sharply with $\mathbf{I}$, which leads to phase-space mixing in $\Omega_{\rm 2D}$. The other frequencies $\omega_{i}$ for $i\neq 1$ are slowly-varying near the saddle point and the distribution along the corresponding angles remain localized over the timescale for phase-space mixing in $\Omega_{\rm 2D}$. In this and the next section we discuss general features of trajectories and observables in the phase space region near a separatrix. We develop this discussion using a simple pendulum, an integrable system with a two-dimensional phase space containing a single saddle point and two separatrices (DLMF, , §22.19).

The Hamiltonian of a simple pendulum is

[TABLE]

where $p$ is the momentum and $\theta\in[-\pi,\pi]$ is the angular position, where $\theta=\pm\pi$ are identical (we have set the pendulum’s length and acceleration due to gravity to one). The point $(\theta,p)=(0,0)$ corresponds to the stable equilibrium, while $(\theta,p)=(\pi,0)$ is its sole saddle point and corresponds to a stationary upright pendulum. Around the saddle point $H_{\rm pend}\sim 2+(p^{2}-x^{2})/2$, where $x=(\theta-\pi)\mod 2\pi$.

Figure 2(a) shows the equal-energy contours in the phase space of the pendulum. Two separatrices, $S+$ and $S-$, divide the phase space into three regions, denoted by $A$, $B$ and $C$, with two distinct kinds of periodic motions: libration and rotation. Libration, confined to region $B$, is an oscillation where $\theta$ is bounded and does not pass the inverted position, $\theta=\pi$. Its time period is $T_{\rm lib}=4K(k)$, where $K(k)$ is the elliptic integral of the first kind DLMF , the modulus $k=\sqrt{\mathcal{E}/2}$ and $\mathcal{E}$ is the energy. Rotation is an unbounded motion in regions $A$ or $C$, where the pendulum passes the inverted position. Its time period is $T_{\rm rot}=2kK(k)$, where $k=\sqrt{2/\mathcal{E}}$. Explicit expressions of libration and rotation motion are given in App. A.

On the separatrices the period is infinite and, hence, the action angle coordinates $(I_{1},\varphi_{1})$ are not defined. Thus, a saddle point precludes the existence of global action-angle coordinates. They are, however, defined separately in each of the three regions. Although, the explicit form of $I_{1}$ and $\varphi_{1}$ in terms of $p$ and $\theta$ is known Brizard (2011), it is not required for our analysis. We will need the location where $\varphi_{1}$ is zero along an equal-energy contour. We define it to be a point near the saddle point where $|p|$ is minimal. This condition is unique for regions $A$ and $C$. In region $B$ there are two such points and we choose the point where $\theta>0$. As the travel time between the two points is a half the period, $\varphi_{1}=\pi$ for the other point. Our choice of $\varphi_{1}=0$ is shown in Fig. 2(a) as dashed-dotted lines originating from the saddle point.

We remark on the properties of solutions on the separatrix, which will be useful later. The two solutions that vary significantly only around $t=0$ and for which $\theta(t=0)=0$ are given by

[TABLE]

Note that $p_{S\pm}(t)$ is well approximated by a bump function (also known as a test function Gelfand and Shilov (1964)) that is nonzero in a finite domain, called the support, and vanishes outside its support. Moreover, an observable $\mathcal{O}(t)$ on the separatrix is (well approximated) by a constant plus a bump function, as long as it is smooth in both $p$ and $\theta$ and periodic in $\theta$.

Trajectories $(\theta(t),p(t))$ that start near one of the separatrices spend most of their time (within a period) near the saddle point as shown with two examples in Fig. 2(b). Changes in $\theta(t)$ and $p(t)$ from their saddle-point value are to good approximation equal to corresponding changes along trajectories on one or more of the separatrices. For example, for the rotation trajectory in Fig. 2(b) the momentum is $p_{A}(t)=p_{\mathrm{S}+}(t-T_{\rm rot}/2)$ for $t\in[0,T_{\rm rot})$, while for the libration trajectory in Fig. 2(b) the momentum is $p_{B}(t)=p_{\mathrm{S}+}(t-T_{\rm lib}/4)+p_{\mathrm{S}-}(t-3T_{\rm lib}/4)$ for $t\in[0,T_{\rm lib})$. In fact, the momentum along any trajectory starting near the saddle point in region $R=A$, $B$ or $C$, respectively, can be written as

[TABLE]

where the sum over $n$ defines the momentum for all $t$ (rather than a single period) and indicator functions $\chi_{R}(s)$ are either zero or one. For the pendulum $\chi_{A}(S+)$, $\chi_{B}(S+)$, $\chi_{B}(S-)$ and $\chi_{C}(S-)$ are one; others are zero. The time shift $t_{0,R}(s)\in[0,T_{R})$ and period $T_{R}$ are determined by the starting point, where $T_{R}=T_{\rm rot}$ and $T_{\rm lib}$ for $R=A,C$ and $R=B$, respectively. Thus, $p_{R}(t)$ is a sum over periodically occurring, non-overlapping bump functions whose support is much smaller than the time period.

The asymptotic symbol $\sim$ in Eq. 10 and elsewhere in this paper implies that either the trajectories start close to the saddle point or the averages are over a Wigner distribution that is initially localized around the saddle-point and whose initial width goes to zero. We also reserve the word asymptotic for these two cases, unless otherwise stated.

VII Long-time expectation value

We now derive the long-time expectation value of observables $\langle\mathcal{O}(t)\rangle$ that are smooth functions of the canonical coordinates $(\mathbf{x},\mathbf{p})$ and depend only on the single action-angle coordinate $\varphi_{1}$ of the subspace $\Omega_{\rm 2D}$ in which the system undergoes phase-space mixing. For periodic coordinates, like angle $\theta$ of the pendulum, we restrict the observables to be periodic in $\theta$. These constraints are not severe as many physically interesting observables have these properties.

The first step is to write the asymptotic form of observable $\mathcal{O}_{R}(t)$ in region $R$, along a trajectory that comes close to the saddle point, in terms of its value along the separatrix trajectories $(x_{1,s}(t),p_{1,s}(t))$ in subspace $\Omega_{\rm 2D}$. Here $s$ labels separatrices. (For the pendulum $s\in\{S+,S-\}$.) We define $\mathcal{O}_{s}(t)=\mathcal{O}(x_{1,s}(t),p_{1,s}(t))$ and realize that $\mathcal{O}_{s}(t)=\mathcal{O}_{\rm sp}+\mathcal{D}_{s}(t)$, where $\mathcal{O}_{\rm sp}$ is the value of the observable at the saddle point and $\mathcal{D}_{s}(t)$ is a bump function localized around $t=0$. Similarly, we decompose $\mathcal{O}_{R}(t)=\mathcal{O}_{\rm sp}+\mathcal{D}_{R}(t)$, where $\mathcal{D}_{R}(t)$ is a series of periodically occurring, non-overlapping bump functions. Then, similar to Eq. 10, we write

[TABLE]

The indicator functions $\chi_{R}(s)$ are system dependent and the sum $s$ is over one or more separatrices.

To compute the long-time limit of $\left\langle\mathcal{O}(t)\right\rangle$ using Eq. 7 we need to evaluate the integral over angle $\varphi_{1}$. (Those over $\varphi_{j}$ for $j>1$ evaluate to unity for allowed observables.) We transform this integral to one over time by choosing a reference trajectory that starts near the saddle point with $\varphi_{1}(0)=0$. For the pendulum, two such trajectories are shown in Fig. 2(b). Then $\varphi_{1}(t)=\omega_{1}t$ and

[TABLE]

For $n=0$ the integrand $\mathcal{D}_{s}(t-t_{0,R}(s))$ is localized around $t=t_{0,R}(s)\in(0,T_{R})$. Its support is enclosed by the integration bounds $t=0$ and $T_{R}$ as the reference trajectory is near the saddle point at these times. For $n\neq 0$ there is no overlap between the support and the integration interval; hence, the integral is zero. We extend the integration limits of $t$ to $(-\infty,\infty)$ for the surviving $n=0$ term and find

[TABLE]

Substituting this expression in Eq. 7 the long-time average becomes

[TABLE]

where the average frequency $\left\langle\omega_{1}\right\rangle_{R}=\int_{R}d\mathbf{I}\,\bar{f}_{R}(\mathbf{I})\,\omega_{1}(\mathbf{I})$ and the expression in the square bracket is independent of the distribution. Equation 14 is an important result of our paper and relates the long-time expectation value of an observable to the mean frequency. The quantity $\mathcal{O}_{\rm sp}$ is the classical value of the observable and the second term is the quantum correction within the TWA.

For the pendulum we assume the initial Gaussian distribution

[TABLE]

where $x=(\theta-\pi)\mod 2\pi$. It is centered around the saddle point, analogous to the Wigner distribution of a mean-field state, where the width $d\ll 1$ 111 The quantum Hamiltonian of a pendulum in the $\theta$ basis is $-(\hbar^{2}/2)\partial^{2}_{\theta}+1-\cos\theta$. The ground state is (approximately) a coherent (Gaussian) state around $\theta=0$ with width $d=\sqrt{\hbar/2}$. When the sign of the potential $\cos\theta$ is suddenly changed, the state becomes dynamically unstable with the initial Wigner distribution as in Eq. 15. . Both $H_{\rm pend}$ and $F_{0}(\theta,p)$ are invariant under the transformations $p\to-p$ and $\theta\to-\theta$. Thus, the time-evolved distribution function is also invariant and observables $\mathcal{O}(\theta,p)$ that are odd functions of either $\theta$ or $p$ have a vanishing expectation value at all times. In contrast, observables that are even functions in both $\theta$ and $p$ can have non-vanishing expectation value.

As an illustration consider $\mathcal{O}(\theta,p)=p^{2}$. Its functional form along the two separatrix solutions in Eq. 9 is the same, i.e., $[p_{\mathrm{S}+}(t)]^{2}=[p_{\mathrm{S}-}(t)]^{2}$ and, using the indicator functions $\chi_{R}(s)$ for the pendulum, we find

[TABLE]

Next we realize that

[TABLE]

where we have used Eq. 9 to evaluate the time integral and defined the “auxiliary frequency” $\varpi$ to be $\omega_{1}$ in region $A$, $C$ and $2\omega_{1}$ in region $B$ with average $\left\langle\varpi\right\rangle=\left\langle\omega_{1}\right\rangle_{A}+\left\langle 2\,\omega_{1}\right\rangle_{B}+\left\langle\omega_{1}\right\rangle_{C}$. From the definition of $\bar{f}_{0,R}(I_{1})$ we also find that

[TABLE]

where the unit-normalized distribution function

[TABLE]

and $\delta(z)$ is the Dirac delta function. The second equality shows that the explicit relationship between $(I_{1},\varphi_{1})$ and $(\theta,p)$ is not required for the analysis.

As shown in App. A.1 the distribution $\mathcal{F}(\varpi)$ is well approximated by a Gaussian when the width $d$ of the initial distribution $F_{0}(\theta,p)$ approaches zero. In fact, the location of its peak value is

[TABLE]

and its width is

[TABLE]

where $\varkappa=0.595\cdots$. Thus, the quantum correction to the long-time expectation value of $p^{2}(t)$ is $1/O(\ln|d|)$.

VIII Time dynamics of relaxation

In this section we study the relaxation of an observable to its long-time expectation value. Observables again depend on only a single angle $\varphi_{1}$ and are periodic in $\varphi_{1}$. We can then write an observable in region $R$ as a Fourier series

[TABLE]

with

[TABLE]

Now, as in Sec. VII, we transform the integral over $\varphi_{1}$ into one over time by choosing a reference trajectory with $\varphi_{1}(0)=0$ and insert $\varphi_{1}(t)=\omega_{1}(\mathbf{I})t$. Using Eq. 11 we find

[TABLE]

where $\delta_{ij}$ is the Kronecker delta, $\alpha_{R}(s)=\omega_{1}t_{0,R}(s)$, the integration variable $\tau=t-t_{0,R}(s)$ and we have suppressed the dependence of $\omega_{1}$ and $T_{R}$ on $\mathbf{I}$. Only the $n=0$ term contributes and

[TABLE]

where the Fourier transform $\mathfrak{D}_{s}\left(x\right)=\int_{-\infty}^{\infty}dt/(2\pi)\,\mathcal{D}_{s}(t)e^{-ixt}$. Substituting this expression into Eq. 22 and using $\mathcal{O}_{R}(t)\equiv\mathcal{O}_{R}(\mathbf{I},\varphi_{1}(t))$ Eq. 5 becomes

[TABLE]

where ${\left\langle\dots\right\rangle}_{R}$ is the average over $f_{0,R}({\mathbf{I}},\bm{\varphi})$, the initial distribution restricted to region $R$. We realize that at long times all Fourier terms except the $m=0$ term must go to zero in order to recover Eq. 14.

We now specialize to the pendulum system. The phases $\alpha_{R}(s)$ are $\alpha_{A}(S+)=\alpha_{C}(S-)=\pi$, $\alpha_{B}(S-)=\pi/2$ and $\alpha_{B}(S+)=3\pi/2$ when $\chi_{R}(s)$ is nonzero and, as shown in App. A.2, we have

[TABLE]

where, as in Sec. VII, the auxiliary frequency $\varpi$ is $\omega_{1}$ in regions $A$, $C$ and $2\omega_{1}$ in region $B$. The distribution $\mathcal{F}(\varpi)$ is well approximated by a Gaussian with mean and width given in Eqs. 20 and 21, respectively. The factor $\varpi\mathfrak{D}_{S+}(m\varpi)$ is slowly varying across the width of $\mathcal{F}(\varpi)$. Carrying out the integral over $\varpi$ in Eq. 27 (after extending the lower limit of the integral to $-\infty$) gives

[TABLE]

Specifically, for $\mathcal{O}(\theta,p)=p^{2}$ we have

[TABLE]

and the time evolution is a sum of oscillatory functions with damping that is Gaussian in time. The oscillation frequency of each term increases linearly with $m$, while simultaneously its damping time, $1/(m\sigma)$, decreases.

IX Condensate in a double-well potential

A Bose-Einstein condensate in a weakly-coupled double-well potential displays Josephson oscillations and macroscopic self-trapping Smerzi et al. (1997); Raghavan et al. (1999); Leggett (2001); Albiez et al. (2005); Zibold et al. (2010). These phenomena are adequately described by a mean-field approximation. Moreover, dynamical instabilities, where quantum effects become important, have also been studied Anglin and Vardi (2001); Chuchem et al. (2010).

A BEC in a symmetric double-well potential is well described by assuming that only two modes $\Psi_{1}(\vec{r})$ and $\Psi_{2}(\vec{r})$ are occupied, one for each well. In the mean-field description the time-dependent order parameter or condensate wavefunction is $\psi_{1}(t)\Psi_{1}(\vec{r})+\psi_{2}(t)\Psi_{2}(\vec{r})$ with complex-valued amplitudes $\psi_{j}(t)$. The real and imaginary parts of $\psi_{j}(t)$ form two pairs of canonical coordinates. Hence, the system has a four-dimensional phase space. Its classical Hamiltonian is

[TABLE]

where $U$ and $J>0$ are the on-site interaction and tunneling energies, respectively Smerzi et al. (1997). The total number of atoms $N=|\psi_{1}|^{2}+|\psi_{2}|^{2}$ and energy $\cal E$ are conserved, making the system is integrable. We note that the underlying quantum Hamiltonian is solvable by the Bethe ansatz Links et al. (2003).

Following the literature it is convenient to introduce $\psi_{j}(t)=\sqrt{N_{j}(t)}e^{i\theta_{j}(t)}$, where $N_{j}$ is the number of atoms in and $\theta_{j}$ is the phase of the condensate in the $j$-th well Smerzi et al. (1997). We can then express Eq. 30 in terms of the fractional population difference $z=(N_{1}-N_{2})/N$ and phase difference $\phi=\theta_{1}-\theta_{2}$, where $\phi\in[0,2\pi]$ and $\phi=0,2\pi$ are identical. In fact, we have $H_{\rm dw}=NJ\times h_{\rm dw}(\phi,z)$, where $h_{\rm dw}(\phi,z)$ is the “single-atom” Hamiltonian that depends on the effective $N$-dependent coupling strength $\Lambda=UN/(2J)$ and is given by

[TABLE]

The Hamiltonian $h_{\rm dw}(\phi,z)$ has a single minimum located at $(\phi,z)=(0,0)$ for $\Lambda>0$. For $\Lambda>1$ the Hamiltonian has a saddle point located at $(\phi,z)=(\pi,0)$. Near the saddle point $h_{\rm dw}(\phi,z)\sim 1+[(\Lambda-1)z^{2}-(\phi-\pi)^{2}]/2$. Figure 3 shows equal-energy contours of $h_{\rm dw}(\phi,z)$ in the two-dimensional phase space $(\phi,z)$ for $\Lambda>1$. Two separatrices $S+$ and $S-$ divide the phase space into regions $A$, $B$ and $C$. Similar to the pendulum, in region $A$ and $C$ the motion is rotational while in region $B$ it is librational. Explicit expressions for rotation and libration trajectories are given in App. B. On each separatrix we consider a trajectory $(\phi_{s}(t),z_{s}(t))$ that only varies significantly around $t=0$ and for which $|z(t)|$ has a maximum at $t=0$. Along these trajectories

[TABLE]

The corresponding $\phi_{S\pm}(t)$ can be calculated by solving $h_{\rm dw}(\phi_{S\pm}(t),z_{\pm}(t))=1$.

We now consider the dynamics of a (zero-temperature) condensate with $N$ atoms prepared at the saddle point within the TWA. We assume that the initial state is $(\psi_{1},\psi_{2})=(\sqrt{N/2},-\sqrt{N/2})$ with corresponding Wigner distribution

[TABLE]

where $i\in\{1,2\}$ and the probability measure is $\prod_{i}d\psi_{i}^{*}\,d\psi_{i}$. The distribution $F_{0}(\psi_{i},\psi^{*}_{i})$ corresponds to the Wigner transform of a product of coherent states, one in each of the two modes with mean atom number $N/2$ and a relative phase of $\pi$.

Observables have a natural interpretation as spin operators when we represent the phase-space $(\phi,z)$ as a sphere with polar angle $\vartheta=\arccos(z)$ and azimuthal angle $\phi$. Hence, observable $z$ corresponds to $s_{z}$, the $z$ component of the unit “spin” $\vec{s}$. The other spin components are $s_{x}=\sin\vartheta\cos\phi=\sqrt{1-z^{2}}\cos\phi$ and $s_{y}=\sin\vartheta\sin\phi=\sqrt{1-z^{2}}\sin\phi$. As in the pendulum case, observables that are odd functions of $\phi$ or $z$ have vanishing expectation values for all times. Thus, $\left\langle s_{z}(t)\right\rangle=\left\langle s_{y}(t)\right\rangle=0$, but $\left\langle s_{x}(t)\right\rangle$ is non-vanishing. Using Eq. 31, we find that $s_{x}=\Lambda z^{2}/2-1$ on the separatrices.

Now we evaluate the long-time limit and time dynamics of $\left\langle s_{x}(t)\right\rangle$. The indicator functions $\chi_{R}(s)$ are $\chi_{A}(S+)=1$, $\chi_{B}(S+)=1$, $\chi_{B}(S-)=1$, $\chi_{C}(S+)=1$ and zero otherwise. Then using Eqs. 14, 32 and following the derivation in Sec. VII we find

[TABLE]

where the auxiliary frequency $\varpi$ is $\omega_{1}$ in regions $A$, $C$ and $2\omega_{1}$ in region $B$. The time evolution of $\left\langle s_{x}(t)\right\rangle$ is found by repeating the steps in Sec. VIII. Details are given in App. B, where we find that the asymptotic expression of $\left\langle s_{x}(t)\right\rangle$ is again given by Eq. 27, with a distribution function $\mathcal{F}(\varpi)$ that is a well approximated by a narrow Gaussian with mean $\mu=\left\langle\varpi\right\rangle$ and width $\sigma\ll\mu$ that depend on $\Lambda$ and $N$. Then Eq. 28 holds and

[TABLE]

It is important to note that, as shown in App. B.1, for large $N$ the mean $\mu$ is $O(1/\ln N)$ and the width is $O[1/(\ln N)^{2}]$. Thus, the quantum correction to the long-time value of $\left\langle s_{x}(t)\right\rangle$ is $O(1/\ln N)$. Quantitative analytical expressions for $\mu$ and $\sigma$ have only been found for $\Lambda-1\ll 1$.

Figures 4(a) and (b) show the long-time expectation value Eq. 34 as a function of $\Lambda$ and Eq. 35 as a function of time, respectively. In addition, the figures show good agreement with numerical TWA results. In the numerical implementation of TWA we sample from the initial distribution $F_{0}(\psi_{i},\psi_{i}^{*})$, propagate the classical equations of motion and compute the expectation value of an observable by averaging over the sample.

X Spinor BEC within the single-mode approximation

A trapped spin-1 (spinor) Bose-Einstein condensate is well-described by a single spatial mode for its three magnetic sublevels Law et al. (1998); Zhou (2001); Zhang et al. (2005). This single-mode approximation (SMA) is valid when the spin healing length, the length scale over which the spin populations of the condensate can change significantly, is larger than the condensate size. The mean-field theory within SMA has turned out to adequately describe atomic spinor experiments with strong spatial confinement Chang et al. (2005); Black et al. (2007); Liu et al. (2009); Bookjans et al. (2011). Quenches to dynamical instability, where quantum effects need to be treated, have also been studied experimentally Hamley et al. (2012); Gerving et al. (2012).

The mean-field wavefunction of the spinor BEC in the SMA is the vector $\vec{\Psi}(\vec{r},t)=\left(\psi_{-1}(t),\psi_{0}(t),\psi_{+1}(t)\right)^{T}\Phi(\vec{r})$, where $\psi_{j}(t)$ is the complex amplitude of the $j$-th magnetic sublevel along the external magnetic field and $\Phi(\vec{r})$ is the time-independent unit-normalized spatial mode. The phase space spanned by the $\psi_{j}(t)$ has six dimensions and the system has three mutually commuting conserved quantities, namely energy, total atom number $N=\sum_{j}|\psi_{j}(t)|^{2}$, and magnetization $M=\sum_{j}j|\psi_{j}(t)|^{2}$. Thus, the system is integrable. We note that the underlying quantum few-mode Hamiltonian is solvable by the Bethe ansatz Bogoliubov (2006); Lamacraft (2011).

It is convenient to write $\psi_{j}(t)=\sqrt{N_{j}(t)}e^{i\theta_{j}(t)}$, where $N_{j}$ and $\theta_{j}$ are the number of atoms in and the condensate phase of sublevel $j$, respectively. Non-trivial dynamics of the spinor system occurs in a reduced two-dimensional space $\Omega_{2D}$ with coordinates $\phi$ and $\rho_{0}$, for a fixed $N$ and $M$. Here, $\phi=\theta_{1}+\theta_{-1}-2\theta_{0}$, where $\phi\in[-\pi,\pi]$ and $\phi=\pm\pi$ are identical; and $\rho_{0}=N_{0}/N$ is the fraction of atoms in the $j=0$ sublevel. In these coordinates the system obeys the “single-particle” classical Hamiltonian Zhang et al. (2005)

[TABLE]

where the coupling strength $c=g_{2}N\int d^{3}r\,|\Phi(\vec{r})|^{4}$ is $N$ dependent, $g_{2}$ is the spin-changing atom-atom interaction strength, the term $q(1-\rho_{0})$ describes atomic level shifts with controllable strength $q$ (in essence due to the quadratic Zeeman interaction) and the conserved unit-magnetization $m=M/N$.

Here, we will only consider a condensate with antiferromagnetic $c>0$ interactions and assume $m=0$. Figure 5 shows equal-energy contours of $h_{\rm spin}(\phi,\rho_{0})$ for a representative $q$ in $(-2c,0)$. The Hamiltonian then has a saddle point at the north pole $\rho_{0}=1$ and $h_{\rm spin}(\phi,\rho_{0})\sim(1-\rho_{0})\{c(1+\cos\phi)+q\}$ with a linear energy dependence for small positive $1-\rho_{0}$. The slope, given in $\{\cdots\}$, changes sign twice when $\phi$ goes from 0 to $2\pi$. Unlike the pendulum and double-well systems, there is only one separatrix $S$, which divides the phase space into regions $A$ and $B$ with rotation and bounded motion, respectively. The expression for $\rho_{0}(t)$ along a general trajectory is given in App. C. The solution along the separatrix that is symmetric about $t=0$ is

[TABLE]

where $y_{1,S}=|q|/(2c)$ and $\Omega=\sqrt{2|q|c(1-y_{1,S})}$. By solving $h_{\rm spin}(\phi_{S}(t),\rho_{0,S}(t))=0$ the corresponding $\phi_{S}(t)$ can be found.

We prepare the system in the mean-field ground state for $q>0$, i.e., $\rho_{0}=1$ or equivalently $(\psi_{+1},\psi_{0},\psi_{-1})=(0,\sqrt{N},0)$. The initial Wigner distribution is

[TABLE]

where $j\in\{+1,0,-1\}$, corresponding to a coherent state for sublevel $j=0$ with a mean atom number $N$ and zero phase and vacuum states for sublevels $j=\pm 1$. The probability measure for the distribution is $\prod_{j}d\psi_{j}^{*}d\psi_{j}$.

The parameter $q$ is then quenched to a value between $-2c$ and [math] at time $t=0$ and the system becomes dynamically unstable. Using Eq. 14 with two contributing regions and one separatrix, the average $\left\langle\rho_{0}(t)\right\rangle$ long after the quench is given by

[TABLE]

where we used the indicator functions $\chi_{A}(S)=\chi_{B}(S)=1$ and defined auxiliary frequency $\varpi$ that is now $\omega_{1}$ in both regions with average $\left\langle\varpi\right\rangle=\left\langle\omega_{1}\right\rangle_{A}+\left\langle\omega_{1}\right\rangle_{B}$. In App. C.1 we show that $\left\langle\varpi\right\rangle\sim 2\pi\Omega/\ln(16N)$. The quantum correction to the long-time value of $\left\langle\rho_{0}(t)\right\rangle$ is, again, $O(1/\ln N)$.

Figure 6(a) shows $\left\langle\rho_{0}(\infty)\right\rangle$ as a function of $q/c$ for $q\in(-2c,0)$ and fixed atom number $N=1000$. The analytical expression of $\left\langle\rho_{0}(\infty)\right\rangle$ with $\varpi=2\pi\Omega/\ln(16N)$ gives a straight line. The figure also shows the predictions from numerical TWA for the same parameters. For small negative $q$ the two curves differ appreciably. We can reproduce the numerical TWA results when we replace $\left\langle\varpi\right\rangle$ in Eq. 39 by its numerical value as obtained from sampling the initial Wigner distribution. For $|q|/c$ much smaller than the scale of our figure, however, the $\left\langle\rho_{0}(\infty)\right\rangle$ from the numerical TWA and that based on computing $\left\langle\varpi\right\rangle$ from sampling still differ. We will return to this issue later on in this section.

The time evolution of $\left\langle\rho_{0}(t)\right\rangle$ is again calculated from Eq. VIII. The dominant contribution to the expectation value is from the trajectories with the action-angle coordinate $\varphi_{1}(0)\approx 0$ (See App. C.2 for a formal justification.) Hence, we can set $\varphi_{1}(0)=0$ and with $\alpha_{A}(S)=\alpha_{B}(S)=\pi$ find

[TABLE]

where $\mathfrak{D}_{S}(x)$ is the Fourier transform of $\rho_{0,S}(t)$. As in the previous section, we define the distribution function $\mathcal{F}(\varpi)$ with $\varpi=\omega_{1}$ in both regions. It is approximately Gaussian with mean $\mu=\left\langle\varpi\right\rangle$ and width $\sigma\ll\mu$ (see App. C). Then, in a manner similar to that used to find Eq. 29, we derive

[TABLE]

Figure 6(b) shows the typical behavior of $\left\langle\rho_{0}(t)\right\rangle$ as a function of time. For long times the evolution is a damped sinusoid oscillating around its asymptotic value, as only one term in the sum significantly contributes. For shorter times the evolution is more complex and multiple terms are important. The numerical TWA simulations are in good agreement with our analytical expression.

At $q=0$ the Hamiltonian $h_{\rm spin}(\phi,\rho_{0})$ has a degenerate line of saddle points along $\phi=\pi$, instead of a single saddle point. The system is then critical and the formalism described so far can not be applied. Nevertheless, we show in App. C.3 that

[TABLE]

where $\alpha=c\,\sqrt{2/N}$ and $F(x)$ is the Dawson integral DLMF . Figure 6(c) shows this evolution as a function of time. The motion seems overdamped with little oscillatory behavior. Agreement with TWA simulation results is very good.

XI Conclusions and outlook

We have analytically studied the time dynamics of two integrable bosonic systems within the truncated Wigner approximation (TWA) when they become dynamically unstable after a quench in a system parameter. The initial Wigner distribution is then centered around a saddle point. We considered a Bose-Einstein condensate (BEC) in a symmetric double-well potential and an antiferromagnetic spinor BEC in the single-mode approximation. Using action-angle variables and the concept of phase-space mixing we derived the long-time expectation value of observables, Eq. 14. We also derived the relaxation dynamics of the expectation value as given in Eq. VIII. We used a simple pendulum as a guide for these derivations.

The time dynamics of the expectation value of an observable is determined by the distribution of frequency $\omega_{1}$ of the classical, periodic trajectories. The evaluation of the time dynamics simplified due to the symmetries of the Hamiltonian and the initial Wigner distribution. These symmetries also motivated the definition of an auxiliary frequency $\varpi$, which has a simple relationship to $\omega_{1}$. For the two bosonic systems when the initial state is a coherent state of $N$ atoms the mean of $\varpi$ is $O(1/\ln N)$. Hence, the deviation of the long-time expectation value from its classical value at the saddle point is $O(1/\ln N)$. The mean determines the typical time scale of the oscillations in the time evolution. The width of $\varpi$ is $O[1/(\ln N)^{2}]$ and determines the relaxation rate. Furthermore, we obtained their explicit dependence on external parameters.

Although we only considered a representative observable for each system, the time dynamics of observables that quantify (condensate) phase or squeezing can be readily computed using our formalism. Our results are also directly applicable to other integrable systems with a single saddle point in phase space, such as a (anti-) ferromagnetic spinor BEC with nonzero magnetization and a BEC in an asymmetric double-well potential. The formalism can be generalized to integrable Hamiltonians with multiple saddle points, for example the Lipkin-Meshkov-Glick model Ribeiro et al. (2008).

We give a brief outlook on the full quantum dynamics of our two bosonic systems and its comparison with the TWA. The Hilbert space of their underlying few-mode quantum Hamiltonians scales linearly with $N$ when restricted to fixed values of conserved quantities. Thus, quenches in these quantum systems can be simulated efficiently on a classical computer. The eigen-energies near the saddle point have been studied using the Wentzel-Kramers-Brillouin (WKB) approximation for a BEC in a double-well potential Nissen and Keeling (2010) (and the Jaynes-Cumming model Berman and Zaslavsky (1995); Keeling (2009); Babelon et al. (2009)). The anharmonicity in the energy-level spacing defines the quantum break time Berman and Zaslavsky (1995), which scales as $O(\ln N)$ near the saddle point Berman and Zaslavsky (1995); Nissen and Keeling (2010). In fact, we find (not discussed here in detail) that the TWA diverges from the quantum dynamics after the first oscillation consistent with this quantum break time. A detailed study will be the subject of a future publication.

Acknowledgements.

This work has been supported by the National Science Foundation Grant No. PHY-1506343. R. M. acknowledges useful discussions with W. Sengupta on phase-space mixing in classical systems.

Appendix A Pendulum

The simple pendulum is used throughout to illustrate our derivation of dynamics and long-time expectation values for few-mode integrable systems. In this appendix we derive results specific to the pendulum. Its Hamiltonian is given in Eq. 8 with canonical coordinates $\theta$ and $p$ satisfying $\{\theta,p\}=1$, where $\{\cdot,\cdot\}$ is the Poisson bracket.

First, librational trajectories $(\theta_{B}(t),p_{B}(t))$ in phase-space region $B$ are DLMF

[TABLE]

where the modulus $k=\sqrt{\mathcal{E}/2}$, $\mathcal{E}$ is the energy of the trajectory and time $t_{0}$ depends on the initial condition. Secondly, rotational trajectories $(\theta_{R}(t),p_{R}(t))$ in regions $R=A$ and $C$ are

[TABLE]

where $k=\sqrt{2/\mathcal{E}}$. The $+$ and $-$ sign correspond to region $A$ and $C$, respectively. The functions $\operatorname{sn}(z,k)$, $\operatorname{cn}(z,k)$ and $\operatorname{dn}(z,k)$ are Jacobi elliptic functions DLMF . Finally, on the separatrices $\mathcal{E}=2$ with trajectories $(\theta_{S\pm}(t),p_{S\pm}(t))$ given by Eq. 9.

A.1 Distribution function $\mathcal{F}(\varpi)$

In this section we calculate the distribution function

[TABLE]

as defined in Eq. 19, as well as its mean and width. Here, the integral is over the whole phase space $\Omega$ and the initial Gaussian distribution $F_{0}(\theta,p)$, given in Eq. 15, has a width $d$ along both $\theta$ and $p$. The auxiliary frequency $\varpi(\theta,p)=\omega_{1}=\pi/[kK(k)]$ in regions $A$ and $C$ and $\varpi(\theta,p)=2\omega_{1}=\pi/K(k)$ in region $B$, where $K(k)$ is the complete elliptic integral of the first kind with modulus $k\in[0,1]$ DLMF .

Near the saddle point the energy $\mathcal{E}\sim 2+(p^{2}-q^{2})/2$, where $q=(\theta-\pi)\mod 2\pi$. The relationship between energy and modulus leads to $k^{2}\sim 1-|p^{2}-q^{2}|/4$ in all regions. Finally, $\varpi\sim\pi/K(k)\sim 2\pi/\ln(64/|p^{2}-q^{2}|)$ using the asymptotic expansion $K(k)\sim\ln(16/k^{\prime 2})/2$ around $k=1$ with complementary modulus $k^{\prime}$ defined by $k^{\prime 2}=1-k^{2}$.

To compute $\cal F(\varpi)$ it is convenient to first introduce the invertible transformation $\mathcal{X}(\varpi)=2k^{\prime 2}/d^{2}\sim 32e^{-2\pi/\varpi}/d^{2}$. The dependence of $\mathcal{X}(\varpi)$ on $d$ will become clear later. We then write

[TABLE]

as $d\to 0$ with the distribution

[TABLE]

and the factor in front of $P(\mathcal{X}(\varpi))$ in the right-hand side of Eq. 48 is the Jacobian $d\mathcal{X}/d\varpi$.

The separatrices divide the neighborhood of the saddle point into four quadrants. We solve Eq. 49 in each quadrant separately. For the quadrant in region $A$ ($p>0$ and $p>q$) we change the integration variables to $p=d\sqrt{2\mathcal{X}}\cosh u$ and $q=d\sqrt{2\mathcal{X}}\sinh u$ with $u\in(-\infty,\infty)$. Similar changes of variables can be used in the other three quadrants (noting that two quadrants lie in region $B$). The contribution to $P(\mathcal{X})$ from each quadrant turns out to be the same and we finally find

[TABLE]

which has no explicit dependence on the width $d$, and $K_{0}(x)$ is a modified Bessel function DLMF . We then have

[TABLE]

as $d\to 0$ and for $\varpi\ll 1$.

Figure 7 shows $\mathcal{F}(\varpi)$ as a function of $\varpi$ for a single $d$. We find that $\mathcal{F}(\varpi)$ is sharply peaked. It approaches zero as $C\,e^{-2\pi/\varpi}/\varpi^{3}$ when $\varpi\to 0^{+}$ and $C$ is a constant. For $\varpi\gtrsim 1$, where Eq. 51 is invalid, either $p$ or $q$ is much greater than $d$ and $F_{0}(\theta,p)$, hence $\mathcal{F}(\varpi)$, is exponentially small. Thus, it is reasonable to approximate $\mathcal{F}(\varpi)$ by a Gaussian as shown in Fig. 7.

We now calculate the mean and variance of $\varpi$ using one of two methods. The mean $\mu\sim\int_{0}^{\infty}d\mathcal{X}\,P(\mathcal{X})\varpi({\mathcal{X}})$ with $\varpi({\mathcal{X}})=2\pi/\ln[32/(d^{2}\mathcal{X})]$. We then identify the small parameter $\lambda^{-1}=\ln(C/d^{2})^{-1}\ll 1$, where the constant $C$ will be determined later, and find

[TABLE]

with the help of the geometric series. Here, $\mathcal{Y}=\ln(C\mathcal{X}/32)$ and ${\rm E}[\mathcal{Y}]$ is the expectation value of $\mathcal{Y}$ with respect to $P(\mathcal{X})$. For $C=64e^{\gamma}$, where $\gamma$ is the Euler-Mascheroni constant, the expectation value $E[\mathcal{Y}]=0$. Hence,

[TABLE]

Similarly, the variance

[TABLE]

and evaluation of the second moment of $\mathcal{Y}$ gives

[TABLE]

Thus, we find $\mu=O(1/|\ln d|)$ and $\sigma=O(1/|\ln d|^{2})$.

The second method estimates $\mu$ and $\sigma$ from the location of and curvature at the maximum of $\mathcal{F}(\varpi)$ using the fact that the distribution is well approximated by a narrow Gaussian. We could not find a closed form for maximum of $\mathcal{F}(\varpi)$. Instead, we present results based on the extremum of $\varpi^{2}\mathcal{F}(\varpi)$. This only introduces small corrections as $\varpi^{2}\mathcal{F}(\varpi)\sim\mu^{2}\mathcal{F}(\varpi)$ over the width of the distribution near $\varpi=\mu$. After some algebra we find

[TABLE]

[TABLE]

where $g(\mathcal{X})=\varpi^{2}\mathcal{F}(\varpi)=4\mathcal{X}K_{0}(\mathcal{X})$ and $\varkappa=0.595\cdots$ is the solution of $dg(\mathcal{X})/d\mathcal{X}=0$.

The estimates of $\mu$ and $\sigma$ obtained by either method gives the same logarithmic scaling with $d$. The numerical prefactors inside the logarithm, however, are different. Figure 7 shows Gaussian distributions with the estimated mean and width based on the two methods. Their difference from the true ${\cal F}(\varpi)$ vanishes as $d\to 0$.

A.2 Time dynamics of observables

In this subsection we derive the time dynamics of observables for a pendulum. That is, we derive Eq. 27 from Eq. VIII. The dependence of the quantity in the bracket $\langle\cdots\rangle_{R}$ in Eq. VIII on the action-angle coordinates is only through $\omega_{1}$ and $\varphi_{1}$. ( This is also true for the other two systems studied in the paper.) Denoting the quantity by ${\cal A}(\omega_{1},\varphi_{1})$ it is then convenient to write

[TABLE]

where

[TABLE]

and $\bm{\varphi}^{\prime}=(\varphi_{2},\dots,\varphi_{n})$ are all the angles except $\varphi_{1}$. (The time dependence of $\cal A$ is suppressed for clarity.) For the pendulum with its 2D phase space Eq. 59 simplifies to $g_{0,R}(\omega_{1},\varphi_{1})=dI_{1}/d\omega_{1}\,f_{0,R}(I_{1},\varphi_{1})$, where $dI_{1}/d\omega_{1}$ is the Jacobian of the transformation between $I_{1}$ and $\omega_{1}$.

The function $g_{0,R}(\omega_{1},\varphi_{1})$ is concentrated around a few points in the $(\omega_{1},\varphi_{1})$ space from the observation that $F_{0}(\theta,p)$ is localized around the saddle point. The justification of this approximation is subtle and technical; it has been relegated to Sec. A.2.1. We find that

[TABLE]

where $\overline{g}_{0,R}(\omega_{1})=\int_{0}^{2\pi}d\varphi_{1}/(2\pi)g_{0,R}(\omega_{1},\varphi_{1})$ is a marginal distribution.

We can now simplify the average and sums on the right-hand side of Eq. VIII into a single average for observables that are even in $\theta$ and $p$. The bump functions $\mathcal{D}_{S+}(t)$ and $\mathcal{D}_{S-}(t)$ are then identical. Moreover, the angular dependence of $g_{0,B}(\omega_{1},\varphi_{1})$ implies that $\left\langle e^{im\varphi_{1}}\right\rangle_{B}=0$ when $m$ is odd so that odd Fourier components in region $B$ do not contribute to $\left\langle\mathcal{O}(t)\right\rangle$. (For regions $A$ and $C$ both even and odd Fourier components contribute.) Using these observations, the definition of the auxiliary frequency $\varpi$ and the values of $\alpha_{R}(s)$, we combine the sum over regions and separatrices into a single sum and arrive at Eq. 27.

A.2.1 Derivation of Eq. 62

We give a quantitative argument for Eq. 62. In the evaluation of $\cal F(\varpi)$ in App. A.1 we observed that each quadrant in the neighborhood of the saddle point contributes equally. In region $A$, where $\varpi=\omega_{1}$, a comparison of Eq. 19 and the definition of $\overline{g}_{0,A}(\omega_{1})$ shows that $\overline{g}_{0,A}(\omega_{1})\propto\mathcal{F}(\omega_{1})$. Thus, $g_{0,A}(\omega_{1},\varphi_{1})$ is localized around $\mu=\left\langle\varpi\right\rangle$ with a width $\sigma\ll\mu$ along the $\omega_{1}$ coordinate.

Next, we define the standard deviation $\Delta_{A}(\omega_{1})$ of $\varphi_{1}$ with respect to the conditional distribution $g_{0,A}(\omega_{1},\varphi_{1})/\overline{g}_{0,A}(\omega_{1})$ at each value of $\omega_{1}$. We now estimate $\Delta_{A}(\omega_{1})$ from the momentum spread $\Delta p_{A}=O(d)$ in region $A$, where $d$ is the width of $F_{0}(\theta,p)$. Using Eq. 46 we find

[TABLE]

where $t_{0}=kK(k)$, because $p_{A}$ is minimal when $\varphi_{1}=0$ (see Fig. 2(a)). Now we expect the relevant $\varphi_{1}$ to be small and use the Taylor expansion $\operatorname{dn}(x+K(k),k)=k^{\prime}+k^{\prime}k^{2}x^{2}/2+\cdots$ for small $x$ to find

[TABLE]

where $k^{\prime}=\sqrt{1-k^{2}}\sim 4e^{-\pi/\omega_{1}}$ and $p_{A}^{\rm min}=2k^{\prime}/k$. Thus, the width $\Delta_{A}(\omega_{1})\propto\omega_{1}\sqrt{\Delta p_{A}/k^{\prime}}\propto\omega_{1}e^{\pi/(2\omega_{1})}\sqrt{d}$. At first glance, this relation contradicts the assumption that $\Delta_{A}(\omega_{1})$ is small because $\Delta_{A}(\omega_{1})$ diverges as $\omega_{1}\to 0^{+}$. From Sec. A.1, however, we know that $\mathcal{F}(\omega_{1})$ and, thus, $\overline{g}_{0,A}(\omega_{1})$ go to zero rapidly as $\omega_{1}\to 0^{+}$. In fact, at the mean value $\omega_{1}=\mu$, given in Eq. 56, we find $\Delta_{A}\left(\mu\right)=O(1/|\ln d|)\ll 1$. Furthermore, $\Delta_{A}(\omega_{1})$ remains small where $g_{0,A}(\omega_{1},\varphi_{1})$ is significant as $\sigma\ll\mu$. Hence, $g_{0,A}(\omega_{1},\varphi_{1})$ is localized in both $\omega_{1}$ and $\varphi_{1}$. (The distribution $f_{0,A}(I_{1},\varphi_{1})$ is not localized in $\varphi_{1}$ as it does not approach zero as $\omega_{1}(I_{1})\to 0^{+}$.)

The nonzero, albeit small, width of $g_{0,R}(I_{1},\varphi_{1})$ in the $\omega_{1}$ coordinate leads to mixing in $\varphi_{1}$. On the other hand, the distribution over $\omega_{1}$ is invariant in time. We can then replace the narrow distribution $g_{0,A}(\omega_{1},\varphi_{1})$ along $\varphi_{1}$ by a delta function. That is, $g_{0,A}(\omega_{1},\varphi_{1})\approx 2\pi\bar{g}_{0,A}(\omega_{1})\delta(\varphi_{1})$. A similar analysis in regions $B$ and $C$ leads to the other two equations in Eq. 62.

Appendix B A condensate in a double-well potential

In this section we derive results pertaining to a two-mode Bose-Einstein condensate in a double-well potential. Its “single-particle” Hamiltonian $h_{\rm dw}(\phi,z)$ is defined in Eq. 31 and $\{\phi,z\}=1$. For $\Lambda>1$ the Hamiltonian has a single saddle point and two separatrices $S$ divide the phase space into three distinct regions $R=A$, $B$, and $C$. The solutions to the equations of motion are Raghavan et al. (1999)

[TABLE]

where

[TABLE]

The “single-particle” energy of the trajectory is $\mathcal{E}$ and $t_{0}$ depends on the initial condition. The corresponding $\phi_{R}(t)$ can be obtained by solving $h(\phi_{R}(t),z_{R}(t))=\mathcal{E}$. (Note that Ref. Raghavan et al. (1999) misses a factor of $1/2$ in the first argument of both $\operatorname{cn}(z,k)$ and $\operatorname{dn}(z,k)$.) Finally, on the separatrices $\mathcal{E}=1$, $\kappa=1$ and $\mathcal{C}=2\sqrt{\Lambda-1}/\Lambda$ with solutions $z_{S\pm}(t)$ given by Eq. 32.

B.1 Distribution function $\mathcal{F}(\varpi)$

We now compute the distribution function $\mathcal{F}(\varpi)$ for a Bose condensate in a double-well potential. The initial Wigner distribution Eq. 33 is localized around the saddle point $(\psi_{1},\psi_{2})=(\sqrt{N/2},-\sqrt{N/2})$. It is convenient to introduce real coordinates $p_{i}$ and $q_{i}$ defined by $p_{1}+iq_{1}=\psi_{1}-\sqrt{N/2}$ and $p_{2}+iq_{2}=\psi_{2}+\sqrt{N/2}$. In these coordinates the Wigner distribution becomes

[TABLE]

where $i\in\{1,2\}$ and the probability measure is $dp_{1}dq_{1}dp_{2}dq_{2}$. Near the saddle point

[TABLE]

and their substitution into $h_{\rm dw}(\phi,z)$ gives the energy

[TABLE]

close to one.

Next, we express the auxiliary frequency $\varpi=\omega_{1}$ in regions $A$, $C$ and $2\omega_{1}$ in region $B$ in terms of coordinates $p_{i}$ and $q_{i}$. From Eq. 65 and the periodicity of elliptic functions, it follows that near the separatrix $\varpi\sim\pi\sqrt{\Lambda-1}/K(k)\sim 2\pi\sqrt{\Lambda-1}/\ln(16/k^{\prime 2})$ where $k=\kappa$ in region $B$ and $1/\kappa$ in regions $A,C$. The modulus $k$ and its complement $k^{\prime}$ depend on $\mathcal{E}$ and thus on the $p_{i}$ and $q_{i}$. With the help of Eqs. 67 and 71, we find

[TABLE]

This choice of $\mathcal{X}$, in particular its $N$ dependence, will simplify later derivations. We realize that $\varpi\sim 2\pi\sqrt{\Lambda-1}/\ln[32N(\Lambda-1)^{2}/(\mathcal{X}\Lambda^{2})]$ and $\mathcal{X}(\varpi)=32N(1-\Lambda^{-1})^{2}e^{-2\pi\sqrt{\Lambda-1}/\varpi}$. Thus, we have established a relation between $\varpi$ and $p_{i}$, $q_{i}$ via the variable $\mathcal{X}$.

The distribution $\mathcal{F}(\varpi)$ is then

[TABLE]

where

[TABLE]

with $\mathcal{Z}(p_{i},q_{i})$ equal to the right-hand side of Eq. 72 and the factor multiplying ${\mathcal{P}}(\mathcal{X})$ in Eq. 73 is the Jacobian $d\mathcal{X}/d\varpi$.

We simplify the integrals in Eq. 74 by changing to “center of mass” and “relative” coordinates $P=(p_{1}+p_{2})/2$, $p=p_{1}-p_{2}$, $Q=(q_{1}+q_{2})/2$ and $q=q_{1}-q_{2}$. We find

[TABLE]

which yields

[TABLE]

Figure 8 shows $\mathcal{F}(\varpi)$ for $N=1000$ and $\Lambda=2$. It is evident from the figure that $\mathcal{F}(\varpi)$ is well approximated by a Gaussian distribution. The mean $\mu$ and width $\sigma$ of $\mathcal{F}(\varpi)$ can be computed from Eqs. 52 and 54, respectively, with $\lambda=\ln[32N(1-\Lambda^{-1})^{2}/\Lambda^{2}]/\sqrt{\Lambda-1}$. Although we have not been able to evaluate analytically the moments $E[\mathcal{Y}^{n}]$ with $\mathcal{Y}=\ln(\mathcal{X})$, the equations imply that $\mu$ is $O(1/\ln N)$ and $\sigma$ is $O[1/(\ln N)^{2}]$.

We can compute $\mu$ using the second method described in Sec. A.1. The location of the maximum of Eq. 73 is a solution to a transcendental equation that does not have a closed form for arbitrary values of $\Lambda$. For small positive $\Lambda-1$, however, we find a closed-form solution by replacing $\cosh$ in Eq. 76 by a constant, chosen such that the approximate $P(\mathcal{X})$ remains unit normalized. Thus,

[TABLE]

and we find

[TABLE]

[TABLE]

where $\varkappa=0.595\cdots$ and $\Lambda-1\ll 1$.

B.2 Time dynamics of observables

The structure of the phase space of a condensate in double-well potential is similar to that of the pendulum. Therefore, we can directly apply the analysis of time dynamics for a pendulum given in Sec. A.2. In particular, the distribution functions $g_{0,R}(\omega_{1},\varphi_{1})$, as defined in Eq. 59, are localized and are given by Eq. 62. Furthermore, observable $\left\langle s_{x}(t)\right\rangle$ obeys Eq. 28.

Appendix C Spinor gas in single-mode approximation

In this section we obtain results for an antiferromagnetic ($c>0$) spinor condensate under SMA. Its “single-particle” Hamiltonian $h_{\rm spin}(\phi,\rho_{0})$ is given in Eq. 36 and $\{\phi,\rho_{0}\}=1$. For $-2c<q<0$ the Hamiltonian has a single saddle point and a separatrix $S$ dividing the phase space into regions $R=A$ and $B$. In both regions Zhang et al. (2005)

[TABLE]

where $\operatorname{cn}(z,k)$ is a Jacobi elliptic function DLMF and $y_{1}\leq y_{2}\leq y_{3}$ are the three real roots of the cubic equation in $\rho_{0}$

[TABLE]

Here, $\mathcal{E}$ is the “single-particle” energy of the trajectory and $m$ is the unit magnetization. In terms of these roots, $\Omega=\sqrt{2|q|c(y_{3}-y_{1})}$ and the modulus $k=\sqrt{(y_{2}-y_{1})/(y_{3}-y_{1})}$. The solution is periodic in time with period $T=2K(k)/\Omega$ and frequency $\omega_{1}=2\pi/T=2\pi\Omega/[2K(k)]$. The corresponding $\phi_{R}(t)$ is obtained by solving $h_{\rm spin}(\phi_{R}(t),\rho_{0,R}(t))=\mathcal{E}$.

On the separatrix $S$ the energy $\mathcal{E}=0$ and the roots of Eq. 81 are $y_{1,S}=|q|/(2c)$ and $y_{2,S}=y_{3,S}=1$. Using the fact $\operatorname{cn}(x,k)\sim\operatorname{sech}(x)$ as $k\to 1$ and setting $t_{0}=0$, we find the separatrix solution

[TABLE]

where $\Omega_{S}=\sqrt{2|q|c(1-y_{1,S})}$.

C.1 Distribution function $\mathcal{F}(\varpi)$

We now study the distribution $\mathcal{F}(\varpi)$ for the spinor condensate by relating the auxiliary frequency $\varpi$ to the conserved quantities $\cal E$, $m$ and $N$. As the initial Wigner distribution $F_{0}(\psi_{j},\psi_{j}^{*})$ is localized near the saddle point with $\rho_{0}=1$, i.e., $(\psi_{+1},\psi_{0},\psi_{-1})=(0,\sqrt{N},0)$, we again define real coordinates $p_{j}$ and $q_{j}$ via $\psi_{j}=\delta_{j0}\sqrt{N}+p_{j}+iq_{j}$. Then the relevant trajectories have energy $\mathcal{E}=0+\widetilde{\mathcal{E}}/N+O(N^{-3/2})$ and unit-magnetization $m=0+{\widetilde{m}}/N+O(N^{-3/2})$, both close to zero. The quantities $\widetilde{\mathcal{E}}$ and ${\widetilde{m}}$ are $O(1)$ and depend on $p_{j}$ and $q_{j}$. We solve for the roots $y_{i}$ perturbatively with small parameter $1/N$ and find that the modulus $k$ is close to one. Then the auxiliary frequency $\varpi=\omega_{1}\sim 2\pi\Omega/\ln(16/k^{\prime 2})$ in regions $A$ and $B$. We define

[TABLE]

which is independent of $N$, and $\alpha=2|q|(1-y_{1,S})/c$. Conversely, $\varpi=2\pi\Omega/[\ln(16N/\mathcal{X})]$. Unlike for the previous two systems, we have not been able to find an analytical expression for the distribution of $\mathcal{X}$. Nevertheless, we can apply Eq. 52 with small parameter $\lambda^{-1}=\Omega/\ln(16N)$ and find

[TABLE]

Moreover, Eq. 54 implies that $\sigma=O[1/(\ln N)^{2}]$; hence, $\sigma\ll\mu$ as $N\to\infty$.

We have numerically evaluated $\mathcal{F}(\varpi)$ and found that it is a Gaussian to a good approximation for $-2c<q<0$. Figure 9 shows $\mathcal{F}(\varpi)$ for $q/c=-1$ and $N=1000$ and a Gaussian fit to this distribution. For fig. 6 we use the mean and width of the numerically obtained $\mathcal{F}(\varpi)$.

C.2 Time dynamics of observables for $-2c<q<0$

We now obtain an approximation for $g_{0,R}(\omega_{1},\varphi_{1})$, as defined in Eq. 59, for the spinor system, where $R\in\{A,B\}$. The initial Wigner distribution $F_{0}(\psi_{i},\psi_{i}^{*})$ is localized around the saddle point and, thus, we expect $g_{0,R}(\omega_{1},\varphi_{1})$ to be localized around the $\varphi_{1}=0$ (see Fig 5). This can be formally justified by writing $\rho_{0}(t)$ along a trajectory near the separatrix in terms of the angle $\varphi_{1}$. Then, similar to Sec. A.2.1 we can show that the spread in $\varphi_{1}$ is much smaller than one where $g_{0,R}(\omega_{1},\varphi_{1})$ is significant. Thus,

[TABLE]

where $\overline{g}_{0,R}(\omega_{1})=\int_{0}^{2\pi}d\varphi_{1}/(2\pi)g_{0,R}(\omega_{1},\varphi_{1})$ is a marginal distribution.

C.3 Time dynamics for $q=0$

The dynamics of a spinor condensate quenched to $q=0$ is qualitatively different from that for $q<0$. Instead of a single saddle point, the Hamiltonian has a degenerate line of saddle points along $\phi=\pi$. Along a trajectory close to this line $\rho_{0}(t)$ is a sinusoid given by

[TABLE]

where energy $\mathcal{E}\equiv h_{\rm spin}(\phi,\rho_{0})>0$ and $t_{0}$ is determined by the initial condition. This trajectory does not spend a significant fraction of its time period near $\rho_{0}=1$ that violates one of the assumptions under which Eq. 14 was derived.

We can, nevertheless, find an analytical expression for $\left\langle\rho_{0}(t)\right\rangle$ by evaluating the expectation value directly from Eq. 1. The initial Wigner distribution, Eq. 38, is localized around $\rho_{0}=1$, and thus time $t_{0}\approx 0$ for the relevant trajectories. Hence, we only require the distribution function

[TABLE]

Now $\rho_{0}=1$ corresponds to the mean-field state $(\psi_{+1},\psi_{0},\psi_{-1})=(0,\sqrt{N},0)$ and near $\rho_{0}=1$ the Hamiltonian $h_{\rm spin}(\phi,\rho_{0})=c[(p_{+1}+p_{-1})^{2}+(q_{+1}-q_{-1})^{2}]/N+O(N^{-3/2})$, with quadratures $p_{j}$ and $q_{j}$ defined by $\psi_{+1}=p_{+1}+iq_{+1}$ and $\psi_{-1}=p_{-1}+iq_{-1}$. Substituting the Wigner distribution into Eq. 87 and computing the integrals, we find $P(\mathcal{E})\sim Nc^{-1}e^{-N\mathcal{E}/c}$. Finally, averaging Eq. 86 over this distribution yields

[TABLE]

where $\alpha=c\,\sqrt{2/N}$ and $F(x)$ is the Dawson integral DLMF .