A semiclassical theory of phase-space dynamics of interacting bosons

We study the phase-space representation of dynamics of bosons in the semiclassical regime where the occupation number of the modes is large. To this end, we employ the van Vleck-Gutzwiller propagator to obtain an approximation for the Green’s function of the Wigner distribution. The semiclassical analysis incorporates interference of classical paths and reduces to the truncated Wigner approximation (TWA) when the interference is ignored. Furthermore, we identify the Ehrenfest time after which the TWA fails. As a case study, we consider a single-mode quantum nonlinear oscillator, which displays collapse and revival of observables. We analytically show that the interference of classical paths leads to revivals, an effect that is not reproduced by the TWA or a perturbative analysis.


Introduction
The crucial difference between quantum mechanics and a statistical theory based on classical mechanics is the method of computing the transition probability between an initial and a final state [1]. In the classical theory, the transition probability is the sum over probabilities of the paths connecting the two states. In contrast, in quantum mechanics, the transition probability is obtained by first summing the amplitudes of all the connecting paths and then squaring the sum. This procedure leads to interference, a feature absent in the classical theory. An archetypal example of interference is a double-slit experiment in which a beam of particles after passing through two slits forms an oscillating intensity pattern on a screen [1].
The aforementioned difference between the theories can be systematically studied in the semiclassical regime (where the typical action ?ÿ, the reduced Planck's constant). In this regime, a probability amplitude can be approximated by the contributions from a subset of all connecting paths: the classical paths [2,3]. (This is the case with the textbook treatment of the double-slit experiment.) Crucially, within this semiclassical approximation, the transition probability retains interference of paths, albeit classical ones. The role of classical trajectories in quantum dynamics was first elucidated by van Vleck [4]. Later, Gutzwiller extended the van Vleck propagator by including Maslov indices and used it to derive his trace formula [5]. The role of classical paths in quantum mechanics has been extensively studied; for example, in scattering [6], localization [7,8], quantum kicked rotor [9], level statistics [10,11], quantum work [12], the Helium atom [13], quantum transport [14][15][16][17], and quantum revivals [18,19].
In this paper, we study a semiclassical approximation of the phase-space dynamics of interacting bosons in the Wigner-Weyl representation. In this representation, the unitary evolution of an initial quantum state in the Hilbert space is equivalent to the evolution of an initial Wigner distribution in phase space in accordance to the Moyal's equation [20][21][22]. The reduction of the state space from a highdimensional Hilbert space to a lower-dimensional phase space makes the phase-space picture particularly useful for implementing approximations of quantum dynamics. An approximation that is usually made is a mean-field approach. In this case, the distribution is approximated at all times by a delta function whose location is determined by the classical Hamilton's equations. The Gross-Pitaevskii equation and its discrete versions fall under this category.
An improvement over the mean-field description is the truncated Wigner approximation (TWA) [23][24][25], where the initial distribution is extended and is the Wigner transform of a quantum state. The subsequent dynamics of the Wigner distribution is still classical. Equivalently, the Moyal's equation is replaced by the classical Liouville's equation. In the literature, the TWA is sometimes called a semiclassical method even though it lacks interference effects. Quantum corrections to the TWA for interacting bosons were studied by Polkovnikov [26,27] using a perturbation theory with the TWA as its zeroth-order approximation. In particular, a nonlinear oscillator was studied whose quantum dynamics exhibits collapse and revival of coherences. The perturbative analysis describes the initial collapse, with increasing accuracy with the order of the perturbation parameter. It fails to describe revivals in the system because the analysis still lacks interference of classical paths.
We perform a semiclassical analysis of a general Bose system in phase space that incorporates interference of classical paths and makes the comparison with the TWA transparent. In particular, our analysis identifies the Ehrenfest time associated with the TWA as the time when the interference of classical paths becomes important. As a case study, we investigate the nonlinear oscillator and show that the semiclassical dynamics leads to revivals. Recently, others have also applied semiclassical methods to bosons. For example, these methods have been applied to coherent backscattering [28] and autocorrelation functions [29] in the Bose-Hubbard model. In addition, the semiclassical Herman-Kluck propagator has been used to study boson dynamics [30,31].
The remainder of the paper is organized as follows. First, we define the phase space of a bosonic system and the Green's function of a Wigner distribution in sections 2 and 3, respectively. A semiclassical approximation of this Green's function is obtained in section 4 and compared with other approximations in the literature [32,33]. In section 4.1, we find that our semiclassical formalism reduces to the TWA when the interference terms are ignored. Next, we discuss Ehrenfest times associated with the TWA and semiclassical approximation in section 4.2. Subsequently, we apply our formalism to analytically study of a nonlinear oscillator in section 5 and conclude in section 6.
2. Phase-space formulation of a bosonic system A bosonic system with a finite number of modes can be described in terms of annihilation and creation operators a ĵ and a ĵ † , respectively, with = j d 1, ..., , where d is the number of modes. For example, the modes could be the sites of a Bose-Hubbard model or spin components of a singlemode Bose-Einstein condensate. The operators satisfy the commutation relations d = a a , , where δ jk is the Kronecker delta function. To construct the phase space, we first define the quadrature operators = is a Dirac delta function and the integrals are over  d . We construct a phase space by imposing {x i , p j }=δ ij , where .,.
{ } is the Poisson bracket. We will refer to = r x p , ( ) as a phase-space point. In other words, the eigenvalues of the quadrature operators are the cartesian coordinates of the phase space. Note that, the number and phase of a bosonic mode are the polar coordinates in this phase space. Thus, by introducing quadrature operators, we have mapped the kinematics of a many-body boson system with d modes to that of a single particle in d-dimensional position or configuration space.
The Wigner transform [22,34] maps an operator  , a function of a ĵ and a ĵ † or x ĵ and p ĵ , to its Weyl symbol  in the phase space. In fact ) · is the dot product between p and q, and the integral is over the configuration space Î  q d . In particular, the Wigner distribution W t r, ( ) at time t is the Weyl symbol of the density operator r t ( ), up to a factor of p . In the Schrödinger picture, the evolution of the Wigner distribution is governed by the quantum Liouville equation [22]. The expectation value of an operator  at a time t is where the integrals are over the phase space  d 2 . Equivalently, in the Heisenberg picture where W r 0 ( ) is the initial Wigner distribution.

Green's function of the Wigner distribution
The Green's function G t r r , , f i ( )of the Wigner distribution in the Schrödinger picture is defined by [21,33,[35][36][37] , , In a seminal paper on quantum dynamics in phase space, Moyal called G t r r , , f i ( )the 'temporal transformation function' [21]. He derived an expression for G t r r , , f i ( ) in terms of Feynman propagators. We give a short and direct derivation.
The time evolution of the density operator is , where U t ( ) and r 0 are the unitary time-evolution and the initial density operator, respectively. We insert ò ñá = dy y y 1 is the Feynman propagator in the configuration space. Next, we express the initial condition r á ñ y y 1 0 2 |ˆ| in terms of the initial Wigner distribution. To this end, we multiply equation (2), evaluated at t=0 and = r r i , by ¢ e ip q i · and integrate over p i to find We substitute this expression in equation (6) and identify = + ¢ y x q Thus, the exact Green's function of the Wigner distribution involves the product of two Feynman propagators in configuration space. We expect that this product will have interference terms.

Semiclassical approximation of the Green's function
A semiclassical approximation of the Green's function G t r r , , f i ( )was obtained using the Weyl transform of the van Vleck-Gutzwiller propagator by [32,33]. This approach incorporates the 'quantum spread' of a classical path, which appears as sub-Planckian oscillations in the Green's function (whose phase space area < ) about classical paths. Here, we derive a coarse-grained approximation that ignores this quantum spreading but, importantly, keeps the interference of the classical paths. Our approximation reduces to the classical propagation, the TWA, for any Hamiltonian when we ignore interference.
A quantum system is said to be in the semiclassical regime when the typical action (in units of ÿ) that appears in the path integral description of the Feynman propagator is much greater than one. For bosonic modes, this regime corresponds to large occupation numbers. In fact, the semiclassical approximation of the propagator, also known as the van Vleck-Gutzwiller propagator, is [2,5,38,39] where the sum is over all classical paths, indexed by b, that start from position x i and reach x f in time t. The action cl ( ) is the position as a function of time τ of the bth classical path with The number of classical paths contributing to K t x x , , f i SC ( ) can be found by studying the dynamics of the initial Lagrangian manifold  i , the hyperplane = x x i . Classical evolution of each point of  i yields a final manifold  f . Figure 1 shows  i and  f for a two-dimensional example. The final manifold folds at singular positions x f on  f where the number of momenta p f on  f as a function of x f changes. Parts of the manifold  f between these singular regions are 'branches'. The example in figure 1 has the three such branches. Crucially, the final momentum which, in general, is a multivalued function of x f at fixed x i and t, is unique on each branch. Therefore, classical paths connecting x i and x f can be indexed by the branches that intersect the manifold = x x f . It is these branches which contribute to the van Vleck-Gutzwiller propagator in equation (9). For example, for the position x f ). An initial manifold  i , the line = x x i shown in black, evolves into the manifold  f , the curve shown in blue, at time t. For example, points A B C , , , and D on  i are mapped to ¢ ¢ ¢ A B C , , and ¢ D on  f , respectively. The classical path connecting C and ¢ C is also shown. The manifold  f has three branches I, II and III separated by caustics ¢ B and ¢ C . The line = x x f intersects  f thrice; therefore, three paths start on  i and reach position x f at time t.
shown in figure 1 has three paths that contribute to the propagator.
Substitution of van Vleck-Gutzwiller propagator in equation (8) yields the semiclassical approximation to the Green's function .
The expression is cumbersome for our analytical study. To proceed, we assume that in equation (10) only the contributions from small regions Q and ¢ Q around = q 0 and ¢ = q 0, respectively, are important; and, secondly, the Taylor expansion of the action up to linear terms is sufficient in these regions. Here, respectively, are the initial and final momenta of the classical path along which the action is computed (see appendix A for a derivation). We further assume that the extent of the small regions Q and ¢ Q in each direction in position space is much greater than  . (Note that from section 2, both position and momentum have the same units as  .) Furthermore, we approximate Substituting these approximations for S b and  b in equation (10) and interchanging the sum and the integral, we find Implicit in the existence of Q is the assumption that x f is away from the position of any caustics, where two branches meet. The example in figure 1 has two caustics. The integration over q and ¢ q yields functions of p f and ) , whose characteristic widths in momentum space are much less than  . (For 'rectangular' regions Q and ¢ Q we obtain multidimensional sinc functions.) Typically, observables are smooth functions in phase space, i.e. they vary slowly on the scale of  . Moreover, initial states of interest are classical states (coherent states) whose width is of the order of  .
where, for clarity, we suppress the dependence of ; ; , 14 . We now show that the 'diagonal' part of the double sum in equation (13), i.e. when = ¢ b b , is equal to G TWA . To this end, we change the independent variables t x p , , where the sum is over all roots p i b (enumerated by b) of Next, we apply the inverse function theorem, which states that the matrix inverse of a Jacobian is the Jacobian of the inverse mapping, to find . Substituting the expression in equation (15), we arrive at which is the diagonal part of equation (13). Thus, the TWA ignores interference of classical paths [29]. For the special 4 The equation is the multidimensional version of the formula d cases of the harmonic oscillator and free particle, the TWA matches with the quantum motion because only a single path contributes to the sum in equation (13) and, hence, there are no interference terms As a corollary, the TWA is a good approximation up to times where there is only a single trajectory contributing to the Green's function.

Ehrenfest times
An Ehrenfest time is the time scale when an approximation to the quantum motion deviates appreciably from exact evolution [40,41]. In fact, there is a hierarchy of Ehrenfest times based on the approximations to the quantum dynamics [42,43]. For the mean-field approximation, the Ehrenfest time t MF is the time scale when an initially localized Wigner distribution becomes distorted and stretched due to nonlinear (not necessarily chaotic) classical dynamics. From section 4.1, we find that the Ehrenfest time t TWA associated with the TWA occurs when the interference of classical paths becomes important consistent with the observations made in [29]. This time scale is greater than t MF because interference of paths occurs when the Wigner distribution becomes so distorted that it fills up the accessible phase space. From the quantum Liouville equation perspective, t TWA is the time scale when the ÿ-dependent quantum terms lead to significant deviation of the quantum evolution of the Wigner distribution from the classical evolution [44]. Finally, there is t SC , the Ehrenfest time for the breakdown of semiclassical approximation based on van Vleck-Gutzwiller propagator, which is greater than t TWA . Numerical studies have shown that the breakdown occurs when diffraction becomes important [45,46].

Case study: a nonlinear oscillator
We consider a single-mode nonlinear oscillator whose quantum Hamiltonian is The nonlinear oscillator has been studied in experiments with a BEC in a deep optical lattice with negligible tunneling [47] and with photons using Kerr nonlinearity [48]. In these experiments, the initial state is well-described by a coherent state, y a ñ = å ñ The collapse time t col is a few times this time constant, as shown in figure 2, and is much smaller t rev for large N.
In the experiments with a BEC in an optical lattice, threebody effects proportional to a a 3 3 (ˆ) † change the nature of the collapse and revival in an interesting manner [49][50][51]. Semiclassical analysis of bosons in optical lattices with small but finite tunneling has been performed in [29] and the semiclassical results agree well with the complicated quantum collapse and revivals of the autocorrelation function for t  t rev . show á ñ a t | ( ) | and á ñ a t ln| ( ) |, respectively, as a function of time t for an initial coherent state whose mean atom number is four. Exact quantum dynamics (labeled Q) displays collapse and revival of á ñ a t | ( ) | with collapse and revival times t col and t rev , respectively. The mean-field solution, labeled MF, is time independent. The TWA result, equation (24), closely replicates the first collapse but shows no revival and deviates appreciably from the quantum dynamics after a time t TWA . On the other hand, the semiclassical approximation (labeled SC), equation (27), agrees well with the quantum evolution for all times. In panel (a) the semiclassical and quantum curves are indistinguishable.

Dynamics according to the TWA
where ρ 2 =x 2 +p 2 . Hence, classical paths lie on circles in phase space centered around the origin (x, p)=(0, 0). Their oscillation frequencies are ω=Uρ 2 /(2ÿ 2 ). The system is integrable as the phase space is two-dimensional and energy is conserved. The angle of the action-angle coordinates is the polar angle j measured in clockwise direction of motion and evolves as j where j i is the initial angle. Using the definition w = ¶ ¶  I NO , we find that the action coordinate I=ρ 2 /2. Hence, ρ is a constant of motion.
For concreteness, let the initial coherent state y ñ 0 | , with occupation number N 1  , be centered along the p-axis, i.e. a = i N. Its Wigner distribution , the expectation value of â within the TWA. Instead of using the Green's function G TWA of equation (14), it is more convenient to work in the Heisenberg picture. In this picture, the Wigner-Weyl transform of the operator a t ( ) is ( ] and j = 0 along the p-axis. Then using equation (4) and writing W x p , 0 ( ) in polar coordinates, we find ò ò For N 1  , it is sufficient to expand the exponent of W x p , 0 ( ) to second order in ρ and j around the location of the maximum of the Wigner distribution, i.e.
Substituting this expression and w r which matches the initial collapse of the coherent state in equation (20), but has no revival. A comparison of equation (26) with the exact quantum result of equation (19) for the absolute value of á ñ a t ( ) is shown in figure 2. The difference between the initial collapses according to exact quantum dynamics and the TWA is small and is only evident in the log plot of figure 2(b). This difference was calculated using a perturbative analysis in [26]. The perturbative analysis, however, lacks the interference of classical paths; hence, it does not explain the revival at t rev and subsequent collapses and revivals. Figure 2 also shows the mean-field value á ñ a t MF |ˆ( ) |, which is a t | ( )| along the single circular trajectory starting from The classical phenomenon of phase-space mixing explains the collapse of á ñ a t ( ) [52]. For an integrable system, the coarse-grained long-time Wigner distribution is uniformly distributed in the angle coordinates of the action-angle variables. For the nonlinear oscillator, µ j a t e i t ( ) ( ) and its expectation value goes to zero as the Wigner distribution mixes in the angle j. Furthermore, within the TWA, the coarsened Wigner distribution reaches a steady state; hence, there is no revival. This latter observation indicates that quantum interference reverses phase-space mixing and revives the quantum state. The role of interference in quantum revivals has been studied for a particle in quartic and infinitesquare-well potentials in [18,19]. In the next section, we find that applying the semiclassical formalism, indeed, leads to revivals in our bosonic system.

Dynamics according to the semiclassical approximation
The calculation of á ñ a t SĈ ( ) according to the semiclassical approximation is lengthy and has been relegated to appendices B and C. We first calculate the action in terms of the polar angles and winding number of classical paths around the origin in appendix B. We carry out the remainder the calculation in appendix C. Here, we list the main steps: 1. The time evolution of an observable with Weyl symbol  x p , ( ) in the Schrödinger picture is given by equation (3). We replace G r r t , , , in equation (C2) into integrals over the initial and final angles j i and j f , respectively, and a double sum over winding numbers of classical paths around the origin. We also express the observable, G r r t , , f i SC ( ), and the initial Wigner distribution in terms of j i , j f and winding numbers. 3. Next, we note that the classical motion in the phase space is restricted in an annulus of radius N 2 and width of O(1). Then, ) . We make similar approximations for the determinants The initial Wigner distribution, however, varies sharply with ρ and requires a more careful approximation. We then solve the remaining integrals. 5 We use ]. We then replace x p ,ˆby their classical limits to obtain  NO , and ignore the second and third terms in the semiclassical limit N ? 1.
Finally, we find where v is difference of the winding number of the interfering paths. This expression corresponds to a train of localized Gaussians and is invariant under the transformation  + t t p U 2 and v v 1; hence, is periodic with time period p =  t U 2 rev . Figure 2 shows that á ñ a t SĈ ( ) agrees with the exact quantum average y y á ñ t a t ( )|ˆ| ( ) for all times. Interestingly, the Wigner distribution, which is initially localized, delocalizes and spreads over the accessible phase space by collapse time t col . It then intermittently forms regular patterns and returns to the initial distribution at = t t rev . Experiments with photons have observed these patterns [48]. (Instead of the Wigner distribution, the Husimi distribution was measured.) The patterns are also apparent in the expectation values á ñ a t n ( ) for > n 1, which can be calculated using our semiclassical method.
Finally, we discuss the Ehrenfest times of the nonlinear oscillator for initial coherent states. From figure 2(a), we see that the mean-field prediction deviates from the quantum evolution well before the collapse time t col , i.e. t < t MF col . In contrast, the deviation of the TWA from quantum evolution (ignoring exponentially small differences) occurs abruptly after a finite time τ TWA =t rev −t col before the first revival of á ñ a t ( ) . The interference of classical paths starts at t inter when the Wigner distribution fills up the annular accessible phase space, i.e. when paths starting from the initial localized distribution with winding numbers zero and one terminate in the same small region of phase space. We can estimate τ inter by noting that for a coherent state the distribution of classical frequencies ω has a mean  UN and width w D =  U N . Therefore, t p w D 2 inter is of the order of t col and, hence, is much smaller than t TWA . In other words, it takes time for the interference of paths to affect á ñ a t ( ) appreciably. In fact, at t TWA the number of interfering classical paths is of the order of N . On the other hand, the Ehrenfest time τ SC is infinite for the nonlinear oscillator.
The Ehrenfest time t TWA depends on the observable under consideration. For example, for á ñ a t 2 ( ) the collapse and revival times are t rev /2 and t col /2, respectively, and the TWA fails after t t 2 rev col ( ) . Nevertheless, t TWA is still greater than t inter for all observables (that are polynomials in a and a † with a degree smaller than N). We also expect the delay in effects of interference and the dependence of t TWA on the observable to hold true for generic integrable systems (where the dynamics is away from singularities like a saddle point of the classical Hamiltonian). In contrast, in a chaotic system and for motion near a saddle point of an integrable system, the Ehrenfest time t t TWA inter [52][53][54].

Conclusion and outlook
In conclusion, we presented a semiclassical theory of phasespace dynamics of bosons. We derived a semiclassical approximation, equation (10), to the exact Green's function of the Wigner distribution. Crucially, the approximation preserves the quantum interference of classical trajectories. In fact, we have shown that the formalism reduces to the TWA when the interference terms are ignored. Hence, the Ehrenfest time associated with the breakdown of the TWA occurs when interference of classical paths becomes important. As a case study, we examined a single-mode nonlinear oscillator whose exact quantum dynamics exhibits collapse and revival. We investigated the dynamics of an observable of this oscillator using the TWA and our semiclassical formalism. Within TWA, the expectation value of an observable collapses due to phase mixing, and there is no revival. The semiclassical approximation, however, reproduces revivals and accurately matches the exact quantum dynamics for all times. Finally, we comment on the long-time validity of our semiclassical approximation. For the nonlinear oscillator, the semiclassical approach is valid for all times 6 . We expect this to be true for generic integrable systems as they can be quantized by the Einstein-Brillouin-Keller method [55]. The situation, however, is not straightforward for chaotic systems. For example, the semiclassical evolution (based on the van Vleck-Gutzwiller propagator) of an initial wavefunction defined on a Lagrangian manifold, whose Wigner distribution is not localized, breaks down after a time of the order of the Ehrenfest time associated with interference of classical paths [35,56]. For localized initial Wigner distributions, however, numerical studies and heuristic arguments have shown that the van Vleck-Gutzwiller propagator works for rather longer times [45,46,57] and only breaks down due to diffraction. The validity of our semiclassical approach for chaotic systems will require further study.

Appendix A. Derivatives of action
We evaluate the partial derivatives of the action S t x x , , b f i ( ) with respect to the initial and final positions. The action satisfies the Hamilton-Jacobi equation and, in principle, its derivatives are well known [58,59]. Here, we give a derivation for the sake of completeness. For notational simplicity, we assume that the configuration space is one-dimensional; generalization to higher dimensions is straightforward. Consider a classical path t t ) and ends at x p , ). Next, consider another classical path whose position in time, . Then the change in the action is  ) starting from a Wigner distribution initially localized in region Ω 0 shown by the solid blue circle. Polar coordinates r j , ( ), with angle j measured from the p-axis in a clockwise direction, are also shown. The region Ω 0 is centered at r j = N , 2 ,0 ( ) ( ) and has a width of O(1). Trajectories starting from Ω 0 lie within the gray annulus. Two paths with traversal time t and winding number zero that start from x=0 and end at = x x f are shown.
Next, we extend the limits on w and ¢ w to ¥ 0, [ ) and write the sums over w and ¢ w in terms of = + ¢ u w w and = -¢ v w w . We combine the sum over u and the integral over j f by defining j j p p = -+ y u mod 2 f i ( ) , whose range is ¥ 0, [ ). We realize that = Here, we derive equation (C4). We restrict our attention to paths that start from the phase-space region Ω 0 , in which the initial Wigner distribution is concentrated. Figure B1 shows the region Ω 0 for the nonlinear oscillator. The paths starting within Ω 0 lie on the annulus shown in the figure. Now, the winding number of a circular path at a fixed traversal time is a stepwise increasing function of the radius. Let the (timedependent) winding numbers of paths that lie on the inner and outer circles of the annulus be w min and w max , respectively, with  w w min max . For a given winding number, there can be two paths that start from Ω 0 with position x i and reach position x f in time t. Figure B1