{'fraction_non_alphanumeric': 0.027249700782141565, 'fraction_numerical': 0.008489436913741754, 'mean_word_length': 4.681214421252372, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum corrections to the dynamics of interacting bosons: beyond the truncated Wigner approximation

We develop a consistent perturbation theory in quantum fluctuations around the classical evolution of a system of interacting bosons. The zero order approximation gives the classical Gross-Pitaevskii equations. In the next order we recover the truncated Wigner approximation, where the evolution is still classical but the initial conditions are distributed according to the Wigner transform of the initial density matrix. Further corrections can be characterized as quantum scattering events, which appear in the form of a nonlinear response of the observable to an infinitesimal displacement of the field along its classical evolution. At the end of the paper we give a few numerical examples to test the formalism.


I. INTRODUCTION
The huge interest to interacting bosons was stimulated by experimental advances in realization of ultracold Bose-Einstein condensates (BECs) 1,2,3,4,5 . If atoms are a subject to an additional periodic potential coming from optical standing waves then they form perfect bosonic crystals with no dislocations, impurities or other defects. Moreover, the hopping amplitude between the adjacent cites of this crystal can be varied by orders of magnitude simply by tuning the intensity of the laser beam 3,5 . Even the sign and the strength of the interaction can be changed using an external magnetic field 6,7 . So one can experimentally address various phenomena without worrying about complications arising from the unwanted degrees of freedom. For example, in equilibrium, the bosons can undergo a transition from a superfluid to an insulator as the strength of the optical potential is increased 8,9,10,11,12,13 . This transition was directly observed in Ref. [3]. The simplicity of the system and relatively slow decoherence suggest a possibility to use the interacting bosons in quantum computing 14,15 . A very appealing idea is to encode quantum information in different sites of the optical lattice rather than in different internal atomic states 16 . In this way the measurement process and the manipulation of the information seem to be quite straightforward.
Another huge advantage of the cold atomic systems is that one can address both theoretically and experimentally dynamic properties of the interacting atoms. In conventional many-body systems strong and non-adiabatic perturbations generically lead to fast damping due to strong coupling to various bath degrees of freedom. So if one is interested in quantum effects, it is extremely hard to work far from equilibrium. There are no such limitations, however, in the atomic systems. Indeed, in order to obtain extremely low temperatures where the atomic de Broglie wavelength becomes comparable with inter-particle spacing, it is necessary to isolate the system to such extent that it can be really thought of as completely closed. One of the examples of the strongly out of equilibrium behavior occurs when the sign of inter-actions changes from repulsive to attractive, for example using Feshbach resonance 17,18,19 . The important issue for nonequilibrium problems from the theoretical point of view is how to take into account quantum fluctuations and go beyond classical or Gross-Pitaevskii (GP) equations of motion. The latter can not adequately describe dynamical properties near instabilities simply because as the fluctuations grow in time, the initial wavepacket spreads very fast. The obvious generalization of the GP equations is the Bogoliubov's approximation 20 . The latter, together with the perturbative treatment of interactions between quasiparticles, gives a consistent expansion of the partition function for equilibrium problems. On the other hand, in nonequilibrium systems it is not always possible to isolate a unique classical path. In particular, if the classical motion can not be described by the linearized version of the GP equations, the whole Bogoliubov's expansion breaks down. In Ref. [21] we considered a specific problem where neither GP nor Bogoliubov's theory would work. In that paper we studied a condensate, which was adiabatically driven to the point of the instability, and showed that because of quantum fluctuations the system evolves into a macroscopically entangled state. We argued that the correct dynamics can be reproduced within the truncated Wigner approximation (TWA) 22,23,24,25,26,27 . The idea behind TWA is that the evolution is still described by the classical GP equations of motion, but the initial conditions for the classical field are distributed according to the Wigner transform of the initial density matrix. We showed that the TWA naturally arises as the first quantum correction to the classical evolution and it gives asymptotically exact short time behavior for any bosonic system. The main purpose of the present investigation is to show how one can go beyond the truncated Wigner approximation.
There have been some works generalizing GP equations by including the interaction of the condensate with excited bosons which are produced during the collapse of the condensate 28,29 . Another possible alternative to the GP and Bogoliubov's approximations is the conventional Keldysh technique 30 . As any diagrammatic expansion, this one relies on the smallness of interaction both for the classical and the quantum fields 31 . Therefore this tech-nique is not suitable for strongly interacting systems near the classical limit. Also it is not suitable for the shorttime dynamics, since the initial conditions are not explicitly written, but rather absorbed into the quantum propagator in a complicated manner 31 . Quite different class of methods, which came from quantum optics, are based on the solution of the exact Fokker-Plank equation for the density matrix which is written in a coherent state representation 22,23,32 . Because of the over-completeness of the coherent states, such a representation is not unique, neither is the resulting Fokker-Plank equation. For example, in the positive P-representation, the evolution of the density matrix can be mapped into the stochastic classical dynamics 22,23 . Within the Wigner approximation, it is easy to reproduce the correct short-time dynamics and recover TWA, but going further in time becomes a tricky issue. Recently, there have been developed an exact stochastic representation of the evolution equations for the density matrix using the Fock basis, which might be preferable to the coherent state one 33,34 . So far this theory is restricted to a particular class of two-body interactions and it is not completely clear (at least to this author) how to generalize it. Let us also point out that a semiclassical expansion of the dynamics for a system of interacting fermions has been discussed in a number of works by Filinov et. al. 35,36,37 . There the authors suggested a perturbative scheme based on the integral representation of the evolution equations for the density matrix. Describing the dynamics in terms of the latter is certainly a valid concept. However, especially near the classical limit, it might be preferable to study the corrections to the evolution equations themselves. And this is exactly what we are going to do here. Our aim is to develop an expansion similar in spirit to that given in Refs. [35,36,37]. However, we will give a completely different derivation, which is suitable for bosonic systems, and will give a simple intuitive interpretation of the obtained results. The other important difference between the two approaches is that we will not attempt to evolve the density matrix itself, all the dynamics will be ascribed to the observable operators.
Another motivation for studying quantum corrections to the classical equations of motion is the possibility to take into account coupling to other external degrees of freedom, which are usually represented by a thermal bath. Though, as we mentioned above, the notion of a bath is not particularly useful for the condensates because of their isolation from the environment, this coupling certainly provides a strong mechanism of decoherence in most of condensed matter systems. A standard way to add such a coupling into the picture is based on the representation of the bath by a set of harmonic oscillators satisfying certain physical properties. After integrating out the environment degrees of freedom one derives Langevin equations, which are essentially classical equations of motion with extra dissipative and random force terms 38 . Unfortunately, the Langevin equations can be rigorously derived only for noninteracting particles in a harmonic potential linearly coupled to the bath, i.e. in the limit where there are no corrections to the classical equations of motion due to quantum effects. So incorporating the latter into the Caldeira-Leggett picture 38 is another important theoretical challenge. Generalization of the derived evolution equations to open systems will be a subject of our future work.
Before going into the actual calculations let us outline the main results of this paper. The quantum corrections manifest themselves in two ways: (i) the initial conditions for the classical equations of motion are distributed according to the Wigner transform of the initial density matrix; the classical observable is the Weyl symbol of the quantum operator, or equivalently its symmetrized version with fields substituted by c-numbers (see also discussion preceding (18)), (ii) there are quantum scattering processes, which are represented as a nonlinear response of the observable to the infinitesimal transform of the fields along their classical trajectories. We would like to emphasize that (i) alone is equivalent to the truncated Wigner approximation. Amazingly the symmetrized quantum operators and the Wigner transform appear automatically in the path integral approach, where initially all the operators are normal-ordered and there is no obvious reason why the Wigner transform should ever emerge. For some situations it is sufficient to ignore the quantum scattering completely and to consider only the evolution along multiple classical paths. We would like to emphasize that any harmonic action (in particular the Bogoliubov's approximation) is completely described by TWA, i.e. by the classical dynamics with appropriate initial conditions. Moreover, to recover the Bogoliubov's approximation it is sufficient to linearize the classical GP equations around the stationary solution. As we showed in Ref. [21], this linearization is not adequate for problems with unstable dynamics. Quantum scattering modifies the classical trajectories themselves, but this modification is conceptually simple. In the first approximation there is only one scattering event (although distributed over the whole space and the time interval of the evolution), i.e., one has to add a small perturbation only once (for a local in time interaction) during the classical evolution and calculate the (nonlinear) response of the observable at the end. In the second order of perturbation theory there are two scattering events, etc. The whole picture thus resembles the perturbative approach to the ordinary scattering problem in the Feynman path integral representation. Although here we mostly concentrate on the interacting bosons, the results are quite general. They are applicable to any type of evolution, where the classical equations of motion arise as a saddle point of the action. Let us also point out that TWA is equivalent to taking into account all classical vertexes in the conventional Keldysh technique 31 . Each quantum scattering event is equivalent to adding a quantum vertex. The perturbative approach we develop here should be recovered also in the Fokker-Planck master equation for the Wigner function 23 , see also Refs. [35,36,37]. However, to obtain the nonlinear response of the observable, one must consider corrections to the distribution function in terms of functional derivatives with respect to the fields, so that the response is recovered upon integrating by parts the product of the distribution function and the observable. While this must be certainly possible to do, we will see that the functional integral derivation gives a very simple and elegant way of deriving the desired equations.

II. QUANTUM CORRECTIONS TO CLASSICAL DYNAMICS
Let us assume that a system of bosons is described by some hamiltonian: H(a j , a † j , t), which in general depends on time and is expressed as a polynomial in creation and annihilation operators a † and a. It does not matter whether a is a continuous in space or a discrete field sitting on a lattice. For a single band boson Hubbard model, which we keep in mind for specific illustrations, the hamiltonian reads: where J is the tunneling constant between the neighboring sites and U is the on-site interaction strength. Both U and J can explicitly depend on time. We also assume that the initial state of the system is given by a density matrix ρ 0 : where |χ represents some basis and P (χ) is the probability to be in a particular state (this P would coincide with the Glauber P-function 22 if {|χ } are the coherent states). If initially the system is in the pure state, the sum contains only a single term. For the Hamiltonian (1), the classical limit is obtained as the number of bosons per site N tends to infinity. So we will develop the expansion in terms of 1/N , where the latter plays the role of the Plank's constant. To see the analogy explicitly let us define Then, according to standard quantum-mechanical formulas: We had to scale the factor of √ N in (3) to have a unique limit for n j at N → ∞. It is easy to recognize in (4) the conventional commutation relation between, say, the position and the momentum with 1/N playing the role of . The reason why the Plank's constant does not enter (4) explicitly is that the phase φ j does not have a classical analogue, while the phase multiplied by does have one; it can be related to the angular momentum in the rotating systems 20 .
A time dependent expectation value of an arbitrary operator Ω is given by: Here T Kτ denotes time ordering along the Keldysh contour going from τ = 0 to τ = t and then back to τ = 0. The exponent of the operator is understood in the usual sense of an infinite product 39 : where τ q = t q/Q and ∆τ = t/Q. The ordering T K τ requires that the multipliers corresponding to later times are placed closer to the operator Ω. Note that contrary to the derivation of the evolution equations within the Keldysh technique 31 , we ascribe dynamics to the operator Ω rather then to the density matrix ρ. In the coherent state basis (5) reads: Here H(a f,b , a ⋆ f,b , τ ) and Ω(a f,b , a ⋆ f,b , τ ) are the normal ordered hamiltonian H and the observable Ω with operators a and a † substituted by complex numbers a f,b and a ⋆ f,b respectively. The expression above is intentionally written in the discrete form, since we want to take spe-cial care of the boundary effects. Instead of the fields a f and a b propagating forward and backward in time, it is convenient to introduce their classical (ψ) and quantum (η) combinations: a f = ψ + η/2, a b = ψ − η/2. The names "classical" and "quantum" are not accidental since in the classical evolution all trajectories are uniquely defined and the backward path should be exactly identical to the forward one. So a q (t) = 0 comes entirely due to quantum fluctuations. After taking the limit ∆τ → 0 in (7) we derive: where L(ψ, ψ ⋆ , τ ) stands for the classical (GP) differential operator acting on the field ψ(t).
Note again that the operator L as well as the fields ψ and η contain spatial indices which we suppressed in (8) to simplify notations. The products like η ⋆ L in (8) are understood as the appropriate sums: j η ⋆ j L j . Equating (9) to zero gives the classical equations of motion. In particular, if H is given by the normal ordered version of (1) with a j and a † j substituted by ψ j and ψ ⋆ j respectively, then L j [ψ, ψ ⋆ , τ ] = 0 is equivalent to the Gross-Pitaevskii equation: Let us explain in some detail how we arrive from (7) to (8). There are several different contributions to the exponent in (7), which we call S. The first one (S 1 ) comes from the terms, which do not involve the hamiltonian H: Here q denotes the discrete time while the spatial indices are suppressed. The first sum in the continuum limit transforms into the integral: which under the substitutions a f → ψ+η/2, a b → ψ−η/2 and after integrating by parts becomes: In the continuum limit the first and the second terms after the sum in (11) clearly go to zero and the last two read: Combining equations (11) -(14) we derive: The second contribution to the exponent of (7) comes from the terms containing the hamiltonian. However, since all of them are proportional to ∆τ , the continuum limit becomes trivial and does not give extra boundary contributions: And finally the last step in our derivation is the expansion of the expression above in powers of η: We intentionally separated terms linear in η and η ⋆ because they enter the classical equations of motion, while the higher powers contribute quantum corrections. Combining (11), (16), and (17) we recover (8).
Let us discuss how quantum fluctuations enter the classical dynamics. From equation (8) we see there are two boundary and one bulk contributions. Thus, integrating out η(t) results only in the modification of the observable operator: Here the average is taken over η with the measure exp(−|η| 2 /2). Note that the original operator Ω entering (7) must be written in the normal ordered form. It is easy to check that integrating out fluctuations according to (18) is equivalent to rewriting Ω in the symmetrized form and substituting a and a † by ψ and ψ ⋆ . We also note that Ω cl coincides with the Weyl symbol of Ω. Let us give a simple example illustrating this statement choosing Ω to be a number operator: According to (18) we have: We see that Ω cl obtained from the normal-ordered Ω using (18) is equivalent to substituting the symmetrized product of a and a † by ψψ ⋆ in (19). Another boundary contribution originates from the field η 0 corresponding to the initial time. Because of the coupling to ψ 0 , this fluctuation introduces a probability distribution for the classical initial conditions: Note that p(ψ 0 , ψ ⋆ 0 ) is nothing but a Wigner transform of the density matrix in the coherent state representation 22,32 and therefore it is not a positively defined quantity. Sometimes p(ψ 0 , ψ ⋆ 0 ) has a weird nonlocal behavior and the semiclassical limit is achieved in somewhat nonintuitive way 32 , see also discussion given in Ref. [21]. If we ignore corrections to the classical equations of motion coming from the higher powers of the quantum field η or the last exponent in (8), then the time dependence of the observable Ω will be given by: where Ω cl (t, ψ 0 , ψ ⋆ 0 ) is evaluated on the classical field ψ(t) satisfying the initial condition ψ(0) = ψ 0 . This expression constitutes the so called truncated Wigner approximation frequently used in quantum optics 22 .
The deviation of the actual trajectories from the classical paths comes from the bulk quantum fluctuations, which appear in the last multiplier in (8). Those can be interpreted as quantum scattering processes. Clearly they are nonzero only if there is an interaction between the bosons, i.e. H contains terms of the order three or higher in the boson fields a and a † . For simplicity, let us restrict the following analysis only to the case of a two particle short range interaction, which is usually a good approximation for atomic gases (see (1)). After the construction is clear, the generalization on other cases becomes straightforward. So We would like to point out that the unique classical limit at N → ∞ is obtained when is kept to be independent of N 40 . Then the quantum scattering part of the action reads: Because S q contains the third power (or in general also higher powers) of the quantum field η, we can treat it perturbatively so that: where f j (τ ) is an infinitesimal shift of the field ψ j (τ ): (27) Let us briefly outline the derivation of (26). Expanding (8) in powers of η gives the following integrals: where H is the classical Hamiltonian of the nonlinear Schrödinger (or GP) equation, not to be confused with H: Next we use the change of variables: and a simple identity: x m e iαx dx = 2π(−i) m δ (m) (α).
Note that the transform from ψ(τ q ) to f (τ q ) is linear and gives no Jacobian for any interaction. The classical equations of motion are recovered from (30) if f ≡ 0. So nonzero f indeed corresponds to the deviation of trajectories from the classical ones. Clearly each term in the expansion in (26) gives an extra prefactor of 1/N 2 , which is the semiclassical parameter for the boson Hubbard model. To see this, one has to keep in mind that all the fields must be rescaled as ψ j → √ N ψ j 40 , so that their expectation values become of the order of 1 and independent of N . Then each derivative with respect to f or f ⋆ in (26) would bring an extra factor of 1/ √ N . The interpretation of (26) and (33) below becomes transparent from the figure 1. The solid lines there denote classical trajectories and the cross represents a quantum scattering event. The quantum corrections appear as a nonlinear response to the infinitesimal displacement of the field, i.e. they reflect the rigidity of the classical motion. It is straightforward to generalize (26) to the interaction nonlocal in space or time. The only difference is that the field ψ and its derivatives over f and f ⋆ would carry different spatial or time indices. The variables f and f ⋆ in  (26) and (33). Solid lines represent classical trajectories, the cross corresponds to a quantum scattering event.
(26) are treated independently. In numerical evaluations it is more convenient to use: Then (26) becomes: We would like to make a few comments about equations (26) and (33). In our specific example we used the number of bosons per well N as a dimensionless semiclassical parameter. As we discussed above, each term in the expansion in those equations brings an extra factor of 1/N 2 . The absence of anywhere, may be confusing since the classical limit certainly corresponds to → 0. However, this should not be surprising since here and quite often in the atomic physics the Planck's constant is either completely absorbed into energies, which are measured in Hz, or into time. We already argued above that the ultimate reason why does not explicitly appear in our formulas is that the phase does not have a classical counterpart. If we use conventional observables like coordinate and momentum or angle and angular momentum, −1 will appear as a prefactor in the action justifying the saddle-point or classical approximation. In the same way the number of bosons per site N appears as a prefactor in the exponent of (8) after the rescaling ψ → √ N ψ and η → √ N η. So in general any expansion in the powers of η is in fact the expansion in powers of . The other important remark is that at small times the deviation from the classical dynamics due to the quantum scattering behaves as O(f (t)/N 2 ), where f (t) is some function of time, which vanishes at t → 0. This proves that the truncated Wigner approximation, where the quantum scattering is completely ignored gives the exact short-time asymptotical behavior of the full quantum dynamics. This very remarkable result is to be contrasted with those obtained within Keldysh technique, where the short time scales are usually unaccessible.

III. NUMERICAL EXAMPLES
To illustrate the current approach let us consider some specific examples. Since the main purpose of this section is not to address particular problems, but rather numerically test the formalism, we consider a simple case of two coupled condensates where it is straightforward to obtain the exact solution and thus to verify the accuracy of the expansion given in (26) and (33). Formally (33), which we use in practice, is a multidimensional integral which is best evaluated using the Monte-Carlo methods. To find the first quantum correction we apply a small shift to a classical field at a random moment of time for each set of initial conditions, then follow the simultaneous evolution of the original and the shifted trajectories and at the end of the evolution calculate the response of the observable using finite differences. The first example we consider is related to the discussion given in our earlier works 21,40 . Namely, we assume that the two condensates, described by the hamiltonian (1) with j = L, R, were initially uncoupled with their wavefunction being a product of two number states: The sub-indices "L" and "R" correspond to the left and right sites respectively. It can be shown 21 that the Wigner transform of the density matrix corresponding to the state (34) is given by: where L N (x) stands for the Laguerre's polynomial of the order N . Then suddenly at t = 0 the tunneling was turned on and the following evolution of the system is studied. In Refs. [21,40] we showed that the Gross-Pitaveskii and the truncated Wigner approximations are very good for sufficiently short time scales t ≤ t c ∼ N/λ = J/U (we remind again that λ = N J/U is the dimensionless measure of the interaction and the classical limit is achieved at N → ∞ keeping λ =const(N )) and they break down completely for longer times t > t c . As an observable it is convenient to choose a scaled number variance: The classical counterpart of Ω can be found either from direct symmetrization of (36) or using (18): In equations (10) we can always rescale time t → t/J so that the dynamics is completely described by a single dimensionless parameter λ. For the illustration we will choose λ = 1 corresponding to the intermediate interaction strength. In figure 2 we plot the resulting evolution of the number variance for the exact solution, truncated Wigner approximation and the first quantum correction. Although account of a single quantum scattering event does not considerably extend the domain of applicability of the classical description, it allows to determine the time t c , where the TWA breaks down without addressing the exact solution. Similarly, the second correction would show the time scale where the first one breaks down, etc.
Relatively small extension of t c for this particular example should not be surprising. As we showed in Ref. [40], the semiclassical description breaks down, because the phase difference accumulated between the different energy levels becomes comparable to π and the discreteness of the spectrum becomes crucial. On the other hand in the semiclassical description the number and the phase are continuously distributed and it becomes very hard to approximate a discrete sum with essentially chaotic phases by a continuous integral. Another example where the discreteness of the spectrum is crucial is quite opposite to that studied above. Assume that initially the system of bosons was in the noninteracting superfluid state, which is characterized by the product of coherent states: In this example we will not assume that the number of sites is equal to two, to avoid unnecessary complications coming from the conservation of the total number of bosons. Then assume that the tunneling was completely turned off. We will follow the evolution of the expectation value of a in a particular site. This example mimics the experiments by M. Greiner et. al. [4], where collapses and revivals of the condensate were observed in this way. In the absence of tunneling, the time evolution in a single site is described by the hamiltonian: By simple time rescaling we can always choose λ = 1. It is easy to solve the corresponding Schrödinger equation in a number state so that It is also straightforward to obtain analytic results within the truncated Wigner approximation and find the first few quantum corrections. Thus: From the structure of (41 -43) it is clear that the further quantum corrections correspond to the 1/N Taylor expansion of (40) and this indeed can be explicitly verified. However, both the truncated Wigner approximation (a T W A (t)) and the quantum corrections fail to reproduce 2πN periodicity in time of the exact result (41). So we can anticipate that the collapse of the condensate occurring at short time scales can be well described, while the revivals can not and this is indeed the case. The reason is that although the function (40) is analytical in N , one needs to use a number of terms exponentially increasing with time to correctly reproduce a(t) . On the physical level we can argue that the restoration of the coherence comes entirely from the discreteness of the spectrum, so that as in the first example, the continuous semiclassical expansion does not work at long time scales. The expec-tation value of the creation operator a as a function of time for the case of the four bosons per well is plotted in the figure 3. Now let us consider few more examples, where the suggested expansion gives considerably better results. Imagine that initially two non-interacting coupled condensates are a subject to the oscillating in time interaction: The frequency ω = 4 is chosen to be exactly twice that of the Josephson-like oscillations in the noninteracting system, so that we can expect a parametric resonance. This process is thus equivalent to a resonant heating of the condensate. We would like to point out that the Gross-Pitaveskij approximation completely fails to describe such a resonance, because the symmetric state is a classical ground state for any strength of interaction λ ≥ −1 21 . The number conserving Bogoliubov's approximation also would fail, because near the parametric resonance we can expect (see figure 4) that a considerable fraction of particles will be excited. The resulting dependence of the relative number variance for eight bosons per well is plotted in figure 4. Obviously, the truncated Wigner approximation gives a very good description of the evolution and the first quantum correction does even a better job so that it hardly deviates from the exact result.
As the last example, we will choose an adiabatic evolution of the ground state of the two coupled condensates with the increasing interaction: Such a form of λ(t) is chosen to mimic the tunneling amplitude J(t) exponentially decreasing with time, with the time being measured in the units of J (see Ref. [21]).
Here δ plays the role of the adiabaticity parameter. If we again take the number variance as an observable, then it will simply decrease with time. The question which remains however, is whether the truncated Wigner approximation will suffice to correctly describe the evolution or not. In the figure 5 we plot the number variance versus interaction for δ = 0.05 and eight bosons per well. It is obvious that the TWA gives a systematic relative deviation from the exact solution, which grows in time and the first quantum correction considerably improves the agreement.

IV. SUMMARY
Let us now summarize and discuss the derived results. We developed a time-dependent perturbation theory around the classical evolution of the system of interacting bosons. We found the two types of corrections. The first one (which is equivalent to the truncated Wigner approximation) does not affect the evolu-tion equations themselves but requires to take an average over an ensemble of trajectories with initial conditions distributed according to the Wigner transform of the initial density matrix. The observable is the classical counterpart of the symmetrized quantum operator (or its Weyl symbol). Further corrections appear in the form of quantum scattering processes, which manifest themselves as a nonlinear response to the infinitesimal change of the fields along their classical evolution (see (26)). We would like to point out again that the widely used Bogoliubov's approximation, or more generally any expansion of the action up to the second order in the fields is entirely contained in the TWA. However, the latter is not limited to the processes, which involve only stable classical evolution (see also discussion in Ref. [21]).
Although we never discussed here the coupling to the thermal bath (or more generally to the external noise source), it is very straightforward to incorporate this into our picture. For example, in a simple representation of the bath as a set of harmonic oscillators, one obtains a dissipative term, which can be directly added to the GP equations and a quadratic term in the quantum fields, which is equivalent to a random force with a Gaussian distribution 38 . In general, there will be other terms in the action as well, but all of them can be treated perturbatively in the same way as we explained irrespective of their origin. So they will result in some kinds of nonlinear responses to the classical, now stochastic, evolution. Note also, that the effects of the dissipation or random forces will wash out all the quantum scattering processes which occur a long time (longer then the relaxation time) prior to the observation. Clearly the response to an infinitesimal perturbation will decay in time if we have such processes. So the theory can be extended to describe the evolution towards the steady states. One have to be cautious though in the low temperature limit, because generically all the relaxation processes are frozen out and the time scale for the scattering leading to the equilibrium goes to infinity. This implies that as T → 0 more and more quantum scattering events have to be considered to correctly describe the equilibrium. So this theory would give some kind of a high temperature expansion. If one is interested in non-equilibrium dynamics, then the limit of the applicability of the perturbative expansion of the given order will be set either by the time of evolution or by the scattering time in a bath, whatever is shorter. Certainly further careful analysis of this approach for the systems interacting with bath is required and it is a subject of the future work.
Finally let us spend a few words on the numerical implementation of the calculations. The solution of the classical equations themselves is very straightforward. Averaging over the initial conditions can be easily done with Monte-Carlo methods and the convergence time either slowly increases with the size of the system or saturates depending on a particular problem. But in any case it does not grow exponentially in size as would be the case for the full quantum solution. The quantum scattering part requires evaluating a nonlinear response, which might be tricky for a nonlinear system of differential equations, but this part is also straightforward and well controlled numerically. We would like to point out that contrary to exact stochastic schemes, this method starts directly from the classical equations of motion, which may be very complicated themselves and shows how to add quantum corrections step by step. So it certainly must be very efficient if the quantum fluctuations are relatively weak. Thus the description of unstable "cat" dynamics considered in Ref. [21], which is very straightforward in the present scheme can be hardly efficiently achieved using the stochastic equations. We believe this method is a competing alternative both to those based on the conventional diagrammatic technique, which is not very suitable for strongly interacting systems, and to the stochastic methods, which can usually deal with a limited number if degrees of freedom.
The author would like to acknowledge helpful discussions with I. Aleiner, I. Carusotto, A. Clerk, S. Girvin, A. Green, A. Melikidze, S. Sachdev, K. Sengupta. This research was supported by US NSF Grant DMR 0098226.