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

###### 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 , where 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.

## 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 , where and are the canonical position and momentum operators, respectively. Here, we have set and the natural frequency of the oscillator to one. Its mean-field ground state is the phase-space point , where , , and is the average over a quantum state. We make the state dynamically unstable by suddenly changing to the Hamiltonian . Under the mean-field equations of motion a dynamically unstable point is stationary. Thus, and holds for all times. In contrast, quantum evolution under leads to exponential growth in the unstable mode Pethick and Smith (2008). In fact, following the language of quantum optics, leads to single-mode squeezing, where 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 time evolves under classical Hamilton’s equations, in contrast to the mean-field approximation where the evolution of a single phase-space point is studied. Here, and are canonical position and momentum coordinates for a classical mean-field Hamiltonian system with degrees of freedom. The initial distribution, , is the Wigner transform Case (2008) of the prequench quantum ground state or any approximation thereof.

For an observable , we define its evolution along a trajectory with initial conditions . The expectation value of over all trajectories is

(1) |

with measures , and the integral is over all phase space . The distribution satisfies for all in accordance with Liouville’s theorem Arnold (1997).

## Iv Classical integrable systems

In classical mechanics a Hamiltonian system with degrees of freedom is integrable if there exist mutually commuting (with respect to the Poisson bracket) conserved quantities Arnold (1997). Then a trajectory in the dimensional phase-space lies on an -dimensional torus. For an integrable system, the coordinates can be transformed to canonical coordinates called actions and angles , such that Hamiltonian is independent of . Crucially, and correspond to the same phase-space point, where is a vector of integers. In these coordinates, the Hamilton’s equations are

(2) |

for all . The frequencies only depend on . Hence, the actions are conserved quantities and the time evolution of the angles has the simple form

(3) |

where and .

For our Hamiltonian systems action-angle coordinates are not globally defined. Instead, they are defined on disjoint regions of by maps from each such region to , where and are the spaces spanned by the actions and angles, respectively. We then construct distribution functions for with normalization . The latter follows from the fact that canonical transformations have a unit Jacobian. The distribution is periodic in and evolves as , where is the initial distribution. Moreover, Eq. 1 becomes

(4) | |||

(5) |

where is the functional form of the observable in region .

## 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 is integrable, where 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 is a function of and the angle is the polar angle in the plane. Points with different rotate around the origin at different frequencies and the distribution function stretches as shown Fig 1. Eventually, the distribution spreads uniformly and mixes in the compact coordinate , while remaining localized in and .

For a general integrable system, the frequencies depend nontrivially on . Hence, the distribution will eventually mix in . 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 by the time-independent distribution (Lifshitz and Pitaevskii, 1981, §1)

(6) |

Consequently, the expectation value at long times becomes

(7) |

Thus, the long-time expectation value of an observable is given by the average over the accessible phase space weighted by .

## 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 spanned by canonical coordinates and . 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 associated with a trajectory in goes to zero as its starting point approaches the saddle point. In fact, near the saddle point varies sharply with , which leads to phase-space mixing in . The other frequencies for are slowly-varying near the saddle point and the distribution along the corresponding angles remain localized over the timescale for phase-space mixing in . 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

(8) |

where is the momentum and is the angular position, where are identical (we have set the pendulum’s length and acceleration due to gravity to one). The point corresponds to the stable equilibrium, while is its sole saddle point and corresponds to a stationary upright pendulum. Around the saddle point , where .

Figure 2(a) shows the equal-energy contours in the phase space of the pendulum. Two separatrices, and , divide the phase space into three regions, denoted by , and , with two distinct kinds of periodic motions: libration and rotation. Libration, confined to region , is an oscillation where is bounded and does not pass the inverted position, . Its time period is , where is the elliptic integral of the first kind DLMF (), the modulus and is the energy. Rotation is an unbounded motion in regions or , where the pendulum passes the inverted position. Its time period is , where . 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 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 and in terms of and is known Brizard (2011), it is not required for our analysis. We will need the location where is zero along an equal-energy contour. We define it to be a point near the saddle point where is minimal. This condition is unique for regions and . In region there are two such points and we choose the point where . As the travel time between the two points is a half the period, for the other point. Our choice of 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 and for which are given by

(9) |

Note that 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 on the separatrix is (well approximated) by a constant plus a bump function, as long as it is smooth in both and and periodic in .

Trajectories 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 and 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 for , while for the libration trajectory in Fig. 2(b) the momentum is for . In fact, the momentum along any trajectory starting near the saddle point in region , or , respectively, can be written as

(10) |

where the sum over defines the momentum for all (rather than a single period) and indicator functions are either zero or one. For the pendulum , , and are one; others are zero. The time shift and period are determined by the starting point, where and for and , respectively. Thus, is a sum over periodically occurring, non-overlapping bump functions whose support is much smaller than the time period.

The asymptotic symbol 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 that are smooth functions of the canonical coordinates and depend only on the single action-angle coordinate of the subspace in which the system undergoes phase-space mixing. For periodic coordinates, like angle of the pendulum, we restrict the observables to be periodic in . These constraints are not severe as many physically interesting observables have these properties.

The first step is to write the asymptotic form of observable in region , along a trajectory that comes close to the saddle point, in terms of its value along the separatrix trajectories in subspace . Here labels separatrices. (For the pendulum .) We define and realize that , where is the value of the observable at the saddle point and is a bump function localized around . Similarly, we decompose , where is a series of periodically occurring, non-overlapping bump functions. Then, similar to Eq. 10, we write

(11) |

The indicator functions are system dependent and the sum is over one or more separatrices.

To compute the long-time limit of using Eq. 7 we need to evaluate the integral over angle . (Those over for 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 . For the pendulum, two such trajectories are shown in Fig. 2(b). Then and

(12) | ||||

For the integrand is localized around . Its support is enclosed by the integration bounds and as the reference trajectory is near the saddle point at these times. For there is no overlap between the support and the integration interval; hence, the integral is zero. We extend the integration limits of to for the surviving term and find

(13) | ||||

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

(14) |

where the average frequency 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 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

(15) |

where . It is
centered around the saddle point,
analogous to the Wigner distribution of a mean-field state,
where the width
^{1}^{1}1
The quantum Hamiltonian of a pendulum in the basis is . The ground state is
(approximately) a coherent (Gaussian) state around with width . When the sign of the potential is suddenly
changed, the state becomes dynamically unstable with the initial Wigner
distribution as in Eq. 15.
.
Both and are invariant under the
transformations and . Thus,
the time-evolved distribution function is also invariant
and observables that are odd
functions of either or have a vanishing expectation value
at all times. In contrast, observables that are even functions in
both and can have non-vanishing expectation value.

As an illustration consider . Its functional form along the two separatrix solutions in Eq. 9 is the same, i.e., and, using the indicator functions for the pendulum, we find

(16) |

Next we realize that

(17) |

where we have used Eq. 9 to evaluate the time integral and defined the “auxiliary frequency” to be in region , and in region with average . From the definition of we also find that

(18) |

where the unit-normalized distribution function

(19) | |||||

and is the Dirac delta function. The second equality shows that the explicit relationship between and is not required for the analysis.

As shown in App. A.1 the distribution is well approximated by a Gaussian when the width of the initial distribution approaches zero. In fact, the location of its peak value is

(20) |

and its width is

(21) |

where . Thus, the quantum correction to the long-time expectation value of is .

## 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 and are periodic in . We can then write an observable in region as a Fourier series

(22) |

with

(23) |

Now, as in Sec. VII, we transform the integral over into one over time by choosing a reference trajectory with and insert . Using Eq. 11 we find

(24) | ||||

where is the Kronecker delta, , the integration variable and we have suppressed the dependence of and on . Only the term contributes and

where the Fourier transform . Substituting this expression into Eq. 22 and using Eq. 5 becomes

where is the average over , the initial distribution restricted to region . We realize that at long times all Fourier terms except the term must go to zero in order to recover Eq. 14.

We now specialize to the pendulum system. The phases are , and when is nonzero and, as shown in App. A.2, we have

(27) | ||||

where, as in Sec. VII, the auxiliary frequency is in regions , and in region . The distribution is well approximated by a Gaussian with mean and width given in Eqs. 20 and 21, respectively. The factor is slowly varying across the width of . Carrying out the integral over in Eq. 27 (after extending the lower limit of the integral to ) gives

(28) |

Specifically, for we have

(29) |

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 , while simultaneously its damping time, , 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 and are occupied, one for each well. In the mean-field description the time-dependent order parameter or condensate wavefunction is with complex-valued amplitudes . The real and imaginary parts of form two pairs of canonical coordinates. Hence, the system has a four-dimensional phase space. Its classical Hamiltonian is

(30) |

where and are the on-site interaction and tunneling energies, respectively Smerzi et al. (1997). The total number of atoms and energy 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 , where is the number of atoms in and is the phase of the condensate in the -th well Smerzi et al. (1997). We can then express Eq. 30 in terms of the fractional population difference and phase difference , where and are identical. In fact, we have , where is the “single-atom” Hamiltonian that depends on the effective -dependent coupling strength and is given by

(31) |

The Hamiltonian has a single minimum located at for . For the Hamiltonian has a saddle point located at . Near the saddle point . Figure 3 shows equal-energy contours of in the two-dimensional phase space for . Two separatrices and divide the phase space into regions , and . Similar to the pendulum, in region and the motion is rotational while in region it is librational. Explicit expressions for rotation and libration trajectories are given in App. B. On each separatrix we consider a trajectory that only varies significantly around and for which has a maximum at . Along these trajectories

(32) |

The corresponding can be calculated by solving .

We now consider the dynamics of a (zero-temperature) condensate with atoms prepared at the saddle point within the TWA. We assume that the initial state is with corresponding Wigner distribution

(33) |

where and the probability measure is . The distribution corresponds to the Wigner transform of a product of coherent states, one in each of the two modes with mean atom number and a relative phase of .

Observables have a natural interpretation as spin operators when we represent the phase-space as a sphere with polar angle and azimuthal angle . Hence, observable corresponds to , the component of the unit “spin” . The other spin components are and . As in the pendulum case, observables that are odd functions of or have vanishing expectation values for all times. Thus, , but is non-vanishing. Using Eq. 31, we find that on the separatrices.

Now we evaluate the long-time limit and time dynamics of . The indicator functions are , , , and zero otherwise. Then using Eqs. 14, 32 and following the derivation in Sec. VII we find

(34) |

where the auxiliary frequency is in regions , and in region . The time evolution of is found by repeating the steps in Sec. VIII. Details are given in App. B, where we find that the asymptotic expression of is again given by Eq. 27, with a distribution function that is a well approximated by a narrow Gaussian with mean and width that depend on and . Then Eq. 28 holds and

(35) | ||||

It is important to note that, as shown in App. B.1, for large the mean is and the width is . Thus, the quantum correction to the long-time value of is . Quantitative analytical expressions for and have only been found for .

Figures 4(a) and (b) show the long-time expectation value Eq. 34 as a function of 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 , 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 , where is the complex amplitude of the -th magnetic sublevel along the external magnetic field and is the time-independent unit-normalized spatial mode. The phase space spanned by the has six dimensions and the system has three mutually commuting conserved quantities, namely energy, total atom number , and magnetization . 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 , where and are the number of atoms in and the condensate phase of sublevel , respectively. Non-trivial dynamics of the spinor system occurs in a reduced two-dimensional space with coordinates and , for a fixed and . Here, , where and are identical; and is the fraction of atoms in the sublevel. In these coordinates the system obeys the “single-particle” classical Hamiltonian Zhang et al. (2005)

(36) |

where the coupling strength is dependent, is the spin-changing atom-atom interaction strength, the term describes atomic level shifts with controllable strength (in essence due to the quadratic Zeeman interaction) and the conserved unit-magnetization .

Here, we will only consider a condensate with antiferromagnetic interactions and assume . Figure 5 shows equal-energy contours of for a representative in . The Hamiltonian then has a saddle point at the north pole and with a linear energy dependence for small positive . The slope, given in , changes sign twice when goes from 0 to . Unlike the pendulum and double-well systems, there is only one separatrix , which divides the phase space into regions and with rotation and bounded motion, respectively. The expression for along a general trajectory is given in App. C. The solution along the separatrix that is symmetric about is

(37) |

where and . By solving the corresponding can be found.

We prepare the system in the mean-field ground state for , i.e., or equivalently . The initial Wigner distribution is

(38) |

where , corresponding to a coherent state for sublevel with a mean atom number and zero phase and vacuum states for sublevels . The probability measure for the distribution is .

The parameter is then quenched to a value between and at time and the system becomes dynamically unstable. Using Eq. 14 with two contributing regions and one separatrix, the average long after the quench is given by

(39) |

where we used the indicator functions and defined auxiliary frequency that is now in both regions with average . In App. C.1 we show that . The quantum correction to the long-time value of is, again, .

Figure 6(a) shows as a function of for and fixed atom number . The analytical expression of with