Exact results for nonequilibrium dynamics in Wigner phase space
EExact results for nonequilibrium dynamics in Wigner phase space
K. Bencheikh ∗ D´epartement de Physique. Laboratoire de physique quantique et syst`emesdynamiques. Universit´e Ferhat Abbas S´etif-1, Setif 19000, Algeria
L. M. Nieto † Departamento de F´ısica Te´orica, At´omica y ´Optica and IMUVA,Universidad de Valladolid, 47011 Valladolid, Spain (Dated: May 25, 2020)We study time evolution of Wigner function of an initially interacting one-dimensional quantumgas following the switch-off of the interactions. For the scenario where at t = 0 the interactionsare suddenly suppressed, we derive a relationship between the dynamical Wigner function and itsinitial value. A two-particle system initially interacting through two different interactions of Diracdelta type is examined. For a system of particles that is suddenly let to move ballistically (withoutinteractions) in a harmonic trap in d dimensions, and using time evolution of one-body densitymatrix, we derive a relationship between the time dependent Wigner function and its initial value.Using the inverse Wigner transform we obtain, for an initially harmonically trapped noninteractingparticles in d dimensions, the scaling law satisfied by the density matrix at time t after a suddenchange of the trapping frequency. Finally, the effects of interactions are analyzed in the dynamicalWigner function. PACS numbers: 03.75.Ss, 05.30.Fk, 73.21.LaKeywords: Dilute atoms gases, quantum quenching, Wigner functions, contact potentials
I. INTRODUCTION
The description of time evolution properties of quantum many body systems is presently a fundamental topicof research. The amazing development of trapping and cooling techniques have led to experimental realization ofquantum many-body systems consisting of ultracold atomic gases where atoms are confined in traps [1, 2]. Theseartificial many-body systems can be produced and loaded in various geometric traps. Such experimental progressallowed a full control of the external parameters in the Hamiltonian governing the quantum system dynamics. Aninteresting issue in the field of ultracold quantum gases is to study the time evolution of a non-equilibrium situationgenerated through a quantum quench, which consists of a sudden change of the Hamiltonian parameters (for examplea change of the harmonic trap frequency or a change in the interaction strength between the atoms of the gas throughFeshbach resonance). A quantum quench is the easiest way to drive a system to non-equilibrium: the system issupposed to be in its Hamiltonian ground state until time t = 0, when the sudden change of a coupling leads to a newHamiltonian according to which the system evolves for t >
0. On the theoretical side, significant advances have beencarried out in understanding fundamental concepts in the non-equilibrium dynamics of quantum many-body system.Among these ideas is the link between quantum dynamics and quantum chaos [3, 4] and the emergence of a newensemble in Statistical Mechanics called generalized Gibbs ensemble, which is a more general concept than the usualgrand-canonical ensemble and turns out to be a powerful tool in the prediction of relaxation processes for certainintegrable one-dimensional systems [5–11].In the present work we investigate the quantum dynamics of an ultra-cold system of N atoms with equal mass m following an interaction quench from finite to zero interaction strength. We start by writing down the underlyingHamiltonians before and after the interaction quench. For t < H = N (cid:88) i =1 p i m + N (cid:88) i =1 V ( r i ) + 12 N (cid:88) i (cid:54) = j N (cid:88) j =1 v ( r i , r j ) , (1) ∗ Electronic address: [email protected] † Electronic address: [email protected] a r X i v : . [ c ond - m a t . qu a n t - g a s ] M a y where p i = − i (cid:126) ∇ i is the momentum operator of the particle i , V is an external confining potential, and v ( r i , r j ) isa two-body atom-atom interaction. We assume that the system is initially in a quantum many-body state | Φ (cid:105) andits associated initial reduced one-body-density matrix ρ ( r , r (cid:48) ) is defined as (see for instance Ref. [12]) ρ ( r , r (cid:48) ) = N (cid:90) d r d r · · · d r N Φ ( r , r , r , . . . , r N )Φ ∗ ( r (cid:48) , r , r , . . . , r N ) . (2)At t = 0, the interactions are turned-off and the many-body Hamiltonian is given by H = N (cid:88) i =1 p i m + N (cid:88) i =1 V ( r i ) , (3)where we are assuming that the external potential may be different before and after the interaction quench. At times t >
0, the system is in a time dependent quantum state given by | Φ( t ) (cid:105) = e − i (cid:126) Ht | Φ (cid:105) . (4)It should be noted that | Φ (cid:105) is not an eigenstate of the post-quench Hamiltonian H . The latter Hamiltonian describes asystem of independent particles, while | Φ (cid:105) describes the state of an initially interacting particles. Before proceedingfurther, we shall first derive a relationship between the initial ρ ( r , r (cid:48) ), and the time dependent ρ ( r , r (cid:48) ; t ) reducedone-body-density matrices. Using the definition ρ ( r , r (cid:48) ; t ) = N (cid:90) d r d r . . . d r N Φ( r , r , r , . . . , r N ; t )Φ ∗ ( r (cid:48) , r , r , . . . , r N ; t ) , (5)we show in Appendix A the following relation ρ ( r , r (cid:48) ; t ) = (cid:90) (cid:90) U ( r , ξ ; t ) ρ ( ξ , ξ ) U ∗ ( r (cid:48) , ξ ; t ) dξ dξ , (6)where U ( r , ξ ; t ) = (cid:104) r | e − i (cid:126) ( p m + V ) t | ξ (cid:105) is the single particle propagator associated to the post quench Hamiltonian inEq. (3). The above relation describes the time evolution of the one-body density matrix for the considered specificquench of interaction, when the system is suddenly driven from interacting to noninteracting configurations.Let us now come to our main concern, and study the subsequent dynamics of the system in phase space. For thisquench scenario we relate the so-called dynamical Wigner distribution function to its initial value just before thequench. Our interest in the Wigner distribution function [13] is motivated by the fact that it provides a useful tool tostudy various properties of many-body systems. Besides, it is well known that it allows a reformulation of quantummechanics in terms of classical concepts [14–16], and it is also used to generate semi-classical approximations [17, 18].Although other distribution functions exist, the Wigner distribution function has the virtue of its mathematicalsimplicity. Nevertheless, it may take negative values as a manifestation of its quantum nature, and therefore it doesnot represent a true probability but a quasi-probability distribution [17, 19]. Wigner distribution functions have beenused in various contexts, as in cold atomic gases [1, 2], quantum optics [20], quantum information [21], quantum chaos[22], and in the study of non-equilibrium dynamics generated by the perturbation of a Fermi gas system [23]. Also, theWigner distribution function of the noninteracting limit of a Fermi gas system at zero and nonzero temperatures hasbeen the subject of recent studies [24–26]. The Wigner function is defined as the Fourier transform of the one-bodydensity matrix ρ ( r , r ; t ) ≡ ρ ( r + s / , r − s / t ) on the relative coordinate s = r − r , where r = ( r + r ) / d dimensions W ( r , p ; t ) = (cid:90) d s ρ ( r + s / , r − s / t ) e − i p . s / (cid:126) , (7)which is a function of phase space variables ( r , p ) at time t . The inverse transformation reads ρ ( r + s / , r − s / t ) = (cid:90) d p (2 π (cid:126) ) d W ( r , p ; t ) e + i p . s / (cid:126) . (8)Setting s = in Eq. (8), one recovers the spatial local density ρ ( r ; t ) := ρ ( r , r ; t ) = (cid:90) d p (2 π (cid:126) ) d W ( r , p ; t ) , (9)normalized to the total particle number as (cid:82) ρ ( r ; t ) d r = N . Integrating W ( r , p ; t ) over the whole space allows toobtain the momentum distribution n ( p ; t ) = (cid:90) d r (2 π (cid:126) ) d W ( r , p ; t ) , such that (cid:90) n ( p ; t ) d p = N. (10)The structure of the paper is as follows. In Section II we consider the case of an initially untrapped one dimensionalinteracting system subjected to a sudden swich-off of interactions. The resulting ballistic dynamics is examined inphase space. We derive a relationship between the dynamical Wigner phase space density at time t > d = 1), has been the subject of an interesting study [27], where it is suddenly subjected to amodification of the trapping frequency and a relationship between the resulting time dependent Wigner function andits initial value is obtained. In this section we propose a generalization to arbitrary spatial dimensions ( d ≥
1) of thisrelation by using an alternative fully quantum mechanical method which is based on the use of time evolution of theone-body density matrix. As a bonus, by using the inverse Wigner transformation to our generalized relationship inWigner phase space, we derive the scaling law of harmonically trapped noninteracting particles in arbitrary dimension d , between the one-body density matrix at time t and its initial value. For d = 1, our result reduces to the one obtainedin [27], as it should.In order to observe how the dynamical Wigner function following a quench is affected by interactions, we considerphysical systems of particles whose dynamics are governed by scaling laws suffering a non-ballistic expansion, and wecompare the resulting Wigner function with the one obtained for a ballistic expansion. Some final conclusions put anend to the paper in Section IV. II. TIME EVOLUTION OF WIGNER FUNCTION OF AN INITIALLY UNTRAPPED INTERACTINGSYSTEM FOLLOWING A SUDDEN SWICH-OFF OF INTERACTIONS.
In this section let us consider the situation where both external potentials V and V in Eqs. (1) and (3) vanish, andthe study is restricted to one-dimensional interacting particles system. The time evolution of the resulting reducedone-body density matrix after a finite time of free expansion t > ρ ( x , x ; t ) = (cid:90) ∞−∞ (cid:90) ∞−∞ U ( x , ξ ; t ) ρ ( ξ , ξ ) U ∗ ( x , ξ ; t ) dξ dξ . (11)Here U ( x, ξ ; t ) = (cid:104) x | e − itH / (cid:126) | ξ (cid:105) with H = − (cid:126) m d dx . The matrix element U ( x, ξ ; t ) is the free particle Feynmanpropagator in the configuration space, given by [28] U ( x, ξ ; t ) = (cid:114) m πi (cid:126) t exp (cid:16) i m (cid:126) t ( x − ξ ) (cid:17) . (12)If we substitute this expression into Eq. (11) we find ρ ( x , x ; t ) = m π (cid:126) t (cid:90) ∞−∞ (cid:90) ∞−∞ e i m (cid:126) t [ ( x − ξ ) − ( x − ξ ) ] ρ ( ξ , ξ ) dξ dξ , (13)which can be rewritten as ρ ( x , x ; t ) = m π (cid:126) t e i m (cid:126) t ( x − x ) (cid:90) ∞−∞ (cid:90) ∞−∞ e i m (cid:126) t [( ξ − ξ )( ξ + ξ ) − x ξ +2 x ξ )] ρ ( ξ , ξ ) dξ dξ . (14)Now it is more convenient to introduce the center-of-mass and relative coordinates, respectively defined by z =( ξ + ξ ) / z (cid:48) = ξ − ξ . Then we can write ρ ( x , x ; t ) = m e i m (cid:126) t ( x − x ) π (cid:126) t (cid:90) ∞−∞ (cid:90) ∞−∞ e i m (cid:126) t zz (cid:48) e − i m (cid:126) t ( x − x ) z e − i m (cid:126) t ( x + x ) z (cid:48) ρ ( z + z (cid:48) , z − z (cid:48) ) dzdz (cid:48) . (15)In the following we shall examine the time evolution of the above one-body density matrix in the Wigner representation.To proceed with computing the Wigner distribution function, we insert Eq. (15) into (7), and we may write W ( x, p ; t ) = m π (cid:126) t (cid:90) ∞−∞ ds e i m (cid:126) t xs (cid:90) ∞−∞ (cid:90) ∞−∞ e i m (cid:126) t zz (cid:48) e − i m (cid:126) t zs e − i m (cid:126) t z (cid:48) x e − ips/ (cid:126) ρ ( z + z (cid:48) , z − z (cid:48) ) dzdz (cid:48) = m π (cid:126) t (cid:90) ∞−∞ (cid:90) ∞−∞ e i m (cid:126) t zz (cid:48) e − i m (cid:126) t z (cid:48) x ρ ( z + z (cid:48) , z − z (cid:48) ) dzdz (cid:48) (cid:90) ∞−∞ ds e i (cid:126) [ mt x − p − mt z ] s . Carrying out the integration on the variable s , we obtain W ( x, p ; t ) = (cid:90) ∞−∞ (cid:90) ∞−∞ e i m (cid:126) t zz (cid:48) e − i m (cid:126) t xz (cid:48) δ (cid:0) x − ptm − z (cid:1) ρ ( z + z (cid:48) , z − z (cid:48) ) dzdz (cid:48) , (16)which after integration over z reduces to W ( x, p ; t ) = (cid:90) ∞−∞ ρ ( x − ptm + z (cid:48) , x − ptm − z (cid:48) ) e − ipz (cid:48) / (cid:126) dz (cid:48) . (17)But according to Eq. (7), the above integral is nothing but the initial Wigner distribution function at t = 0, and atthe phase space point ( x − pt/m, p ), so that for t > W ( x, p ; t ) = W ( x − ptm , p ) . (18)This relation is the Wigner phase space version of Eq. (11). Upon integration over x and making use of the definitionin Eq. (10), we get for t > n ( p, t ) = (cid:90) ∞−∞ dx π (cid:126) W ( x, p ; t ) = (cid:90) ∞−∞ dx π (cid:126) W ( x − ptm , p ) = (cid:90) ∞−∞ dξ π (cid:126) W ( ξ, p ) = n ( p, . (19)Therefore, we find that under ballistic expansion (where interatomic collisions during the expansion are not present)the dynamical momentum density for t > n ( p, n ( p,
0) exhibits a 1 /p tail at large momentum (see for instance Ref. [29]), equation (19) says that thislong tail behavior is preserved in the dynamical momentum density for all positive times.It is worth noticing that Eq. (11) remains valid for an initial non-vanishing confining potential ( V (cid:54) = 0) providedthat one realizes at t = 0 a simultaneously sudden double quench, where the interactions are turned-off with therelease of the trap ( V = 0). As a consequence the results in Eqs. (18) and (19), in this case hold also true. A. Case of two particles interacting through an attractive δ interaction As a first case study we consider a simple model system of two particles interacting through an attractive Diracdelta potential with highly asymmetric mass imbalance: the particle with mass M is so heavy that the center-of-massmotion can be ignored. The problem reduces then to a one-body system and the m mass light particle Hamiltonianis given for t ≤ H a = − (cid:126) m d dx − aδ ( x ) . (20)Here a > φ ( x ) = √ α e − α | x | with energy E a = − ma / (2 (cid:126) ), and α = ma/ (cid:126) . The corresponding Wigner distribution function for t ≤ W ( x, p ) = 2 α (cid:126) e − α | x | p + α (cid:126) (cid:20) cos (cid:18) p | x | (cid:126) (cid:19) + α (cid:126) p sin (cid:18) p | x | (cid:126) (cid:19)(cid:21) . (21)A plot of this function is given on the left hand side of Figure 1.At time t = 0 the interaction is suddenly turned off ( a = 0) and then, from (18) we get the following Wignerfunction for t > W ( x, p ; t ) = 2 α (cid:126) e − α | x − pm t | p + α (cid:126) (cid:34) cos (cid:32) p (cid:12)(cid:12) x − pm t (cid:12)(cid:12) (cid:126) (cid:33) + α (cid:126) p sin (cid:32) p (cid:12)(cid:12) x − pm t (cid:12)(cid:12) (cid:126) (cid:33)(cid:35) . (22)A plot of this function W ( x, p ; t ) is given for a particular value of t > t = 0 there are two symmetry axes, x = 0 and p = 0; for bigger valuesof t the perfect symmetry is lost, although the axes p = 0 and p = mx/t (this one clockwise rotating with time) playan important role; squeezing is more and more pronounced as t → ∞ , the axis p = mx/t approaching the axis p = 0. FIG. 1: The Wigner function for an attractive Dirac δ interaction obtained in Eq. (22), for the values (cid:126) = m = a = 1. Onthe left for t = 0, where two symmetry axes, x = 0 and p = 0, are evident; on the right after the quenching, for t = 1, a cleardistortion of the initial function can be appreciated: the perfect symmetry is lost and squeezing appears between the axes p = 0and p = x . The white regions correspond to small negative values of W ( x, p ; t ). Finally, using (19) it is easy to show analytically that, the initial momentum density is n ( p,
0) = 2 α (cid:126) π ( p + α (cid:126) ) , (23)It is easy to check that this result satisfies the second equation in (10) for N = 1 particle. It has been shown in [29]that for a system of zero-range δ − interacting one-dimensional atoms with arbitrary strength, the high- p asymptoticbehavior of the momentum distribution for both free and harmonically trapped atoms, exhibits a universal 1 /p dependence. As can be seen in Eq. (23) and at large values of the momentum p , we recovered this 1 /p dependenceof the momentum distribution. B. Case of two particles interacting through a δ - δ (cid:48) interaction We are going to consider now an extension of the previous study of two interacting particles that takes into accountthe presence of an extra point-like interaction term in the potential, proportional to δ (cid:48) . This type of point or zero-range potentials are a subject of recent study in differents contexts [30–33]. The Hamiltonian of the light particlewith mass m is now given for t < H a,b = − (cid:126) m d dx − aδ ( x ) + bδ (cid:48) ( x ) , a > , b ∈ R . (24)The associated Schr¨odinger equation has been carefully analyzed in [30], where it was proven that the above Hamil-tonian supports only one bound state of energy E a,b = − ma (cid:126) (cid:0) m b (cid:126) (cid:1) = − (cid:126) κ m , κ > . (25)The normalized wave function is φ a,b ( x ) = A e − κ | x | (1 + B sign( x )) , (26)where sign( x ) stands for the sign function and A = √ ma (cid:126) (1 + B ) , B = mb (cid:126) . (27)In this case the Wigner distribution function for t ≤ W ( x, p ) = 2 κ (cid:126) e − κ | x | (1 + B )( p + κ (cid:126) ) (cid:20) (1 − B ) cos (cid:18) p | x | (cid:126) (cid:19) + κ (cid:126) + B (2 p + κ (cid:126) ) + 2 B ( p + κ (cid:126) ) sign( x ) κ (cid:126) p sin (cid:18) p | x | (cid:126) (cid:19)(cid:21) . (28)Remark that the presence of the sign function on (28) indicates the presence of a discontinuity of the Wigner function W ( x, p ; 0) along the line x = 0. If we consider the limit B → b = − . b = − . b = − . b = − . FIG. 2: The Wigner function for a Dirac δ - δ (cid:48) interaction as in (24), for the values (cid:126) = m = a = 1. From left to right, and fromtop to bottom: b = − . b = − . b = − . b = − .
2. The white regions correspond to small negative values of W ( x, p )in (28). The discontinuity of the surface at x = 0 can be clearly seen on the plots. At time t = 0 the interaction is turned off ( a = b = 0) and then, from (18) we get the following Wigner functionfor t > W ( x, p ; t ) = 2 κ (cid:126) e − κ (cid:12)(cid:12)(cid:12)(cid:12) x − ptm (cid:12)(cid:12)(cid:12)(cid:12) (1 + B )( p + κ (cid:126) ) (cid:34) (1 − B ) cos (cid:32) p (cid:12)(cid:12) x − ptm (cid:12)(cid:12) (cid:126) (cid:33) (29)+ κ (cid:126) + B (2 p + κ (cid:126) ) + 2 B ( p + κ (cid:126) ) sign( x − ptm ) κ (cid:126) p sin (cid:32) p | x − ptm | (cid:126) (cid:33)(cid:35) . Plots of this function W ( x, p ; t ) are given in Figure 3 for t = 1 and values b = − . b = − . p = 0 is preserved, but the axis x = 0 rotatesaround the origin and becomes x = p (this one changing with time); squeezing is stronger as t → ∞ , the axis x = 0approaching the axis p = 0. FIG. 3: The Wigner function W ( x, p ; t = 1) for a Dirac δ - δ (cid:48) interaction with values (cid:126) = m = a = 1, for b = − . b = − . p = 0 and p = x . The white regions correspond to small negative values of W ( x, p ; t = 1). The discontinuity of the surface at x = p canbe clearly seen. From (19) and (28) the initial momentum density can be determined analytically: n ( p,
0) = 2 (cid:126) απ B (1 + B ) p + (cid:126) α ((1 + B ) p + (cid:126) α ) . (30)Again, it is easy to check that this result satisfies the second equation in (10) for N = 1 particle. In addition itcoincides with (23) in the limit B = mb/ (cid:126) → /p for high p if B (cid:54) = 0 (if B = 0 it was alreadymentioned after (23) that it scales as 1 /p ). Hence, it is interesting to observe that the presence of the additional δ (cid:48) interaction term, changes substantially the high momentum tail of the momentum density: at large values of themomentum p , the momentum distribution in Eq. (30) exhibits 1 /p behavior while in the absence of δ (cid:48) interactionthis density scales as 1 /p . A plot of the momentum density as a function of p and B in Figure 4 clearly shows twoeffects. The aforementioned behavior at large values of p and the appearance of a local minimum around p = 0.This minimum is due to the presence of a polynomial in the numerator of the expression of n ( p ) generated by theadditional δ (cid:48) interaction term. We observe that this minimum is more pronounced as B increases and transforms to amaximum for smaller B values. To be more precise, the point p = 0 is a minimum of n ( p ) for B > B > p = ± (cid:126) α √ B − B (1 + B ) . FIG. 4: Momentum density n ( p ) as given by (30) as a function of p and B ≥
0, with values (cid:126) = m = a = 1. For B = 0 werecover the profile of (23), which is similar for B ≤
2. For values such that B > p = 0: p = 0 turns from a maximum to a minimum and two symmetrical maxima appear at the points p = ± (cid:126) α √ B − / ( B (1 + B )). III. TIME EVOLUTION OF WIGNER DISTRIBUTION FUNCTION FOLLOWING A SUDDENCHANGE OF TRAPPING HARMONIC POTENTIAL IN D DIMENSIONS
Very recently Dean et al. [27] computed the time evolution of the Wigner function W ( x, p ; t ) for a one-dimensionalsystem of noninteracting fermions (called also system of independent fermions), initially trapped by a harmonicoscillator potential with frequency ω , subjected to a sudden change of frequency passing from ω to ω . At time t ,after this quench it was found (see Eq. (13) of [27]) that W ( x, p ; t ) = W (cid:18) x cos ωt − p sin ωtmω , mωx sin ωt + p cos ωt (cid:19) , (31)where W is the initial Wigner function of the system. To derive this relationship the Liouville equation is used. Herewe propose an alternative fully quantum mechanical derivation based on the use of time evolution of one-body densitymatrix which is valid for an arbitrary spatial dimension d , which obviously reduces to the result in (31) for d = 1 asit should.Consider a system of noninteracting fermions initially confined in a d − dimensional isotropic harmonic oscillatorpotential of the form V ( r ) = mω r /
2. Suddenly the frequency is changed from ω to ω and the system expandsin the potential V ( r ) = mω r /
2. According to (6), the time evolution of the resulting one-body density matrix at t > ρ ( r , r ; t ) = (cid:90) (cid:90) U osc ( r , ξ ; t ) ρ ( ξ , ξ ) U ∗ osc ( r , ξ ; t ) dξ dξ , (32)Here ρ ( ξ , ξ ) is the one-body density matrix at t = 0 and U osc ( r , ξ ; t ) is the well known propagator associated tothe isotropic harmonic oscillator potential with frequency ω , given in d dimensions by [28] U osc ( r , ξ ; t ) = (cid:16) mω πi (cid:126) sin ωt (cid:17) d exp (cid:104) i mω (cid:126) sin ωt [( r + ξ ) cos ωt − r · ξ ] (cid:105) (33)by substitution we obtain ρ ( r , r ; t ) = (cid:18) mω π (cid:126) | sin ωt | (cid:19) d (cid:90) (cid:90) ρ ( ξ , ξ ) e imω (cid:126) sin ωt [( r + ξ − r − ξ ) cos ωt − r .ξ +2 r · ξ ] dξ dξ . Using the centre of mass u = ( ξ + ξ ) / v = ( ξ − ξ ) coordinates we can write ρ ( r , r ; t ) = (cid:18) mω π (cid:126) | sin ωt | (cid:19) d e i mω cos ωt (cid:126) sin ωt ( r − r ) (cid:90) (cid:90) ρ (cid:16) u + v , u − v (cid:17) e imω (cid:126) sin ωt [(2 u cos ωt − ( r + r )) . v − u . ( r − r )] d u d v . According to the definition in Eq. (7), the Wigner function of this density matrix can be written as W ( r , p ; t ) = (cid:18) mω π (cid:126) | sin ωt | (cid:19) d (cid:90) (cid:90) e imω (cid:126) sin ωt ( u cos( ωt ) − r ) v ρ (cid:16) u + v , u − v (cid:17) d u d v (cid:90) e i (cid:126) ( mω r cot g ( ωt ) − mω sin ωt u − p ) . s d s , (34)with r = ( r + r ) / s = r − r . We can perform the d dimensional s integration, (cid:90) e i (cid:126) ( mω r cot g ( ωt ) − mω sin ωt u − p ) . s d s = (2 π ) d δ (cid:16) mω (cid:126) r cot g ( ωt ) − mω (cid:126) sin ωt u − p (cid:126) (cid:17) = (2 π (cid:126) ) d (cid:18) | sin ωt | mω (cid:19) d δ (cid:18) r cos ωt − u − sin ωtmω p (cid:19) (35)where we have used the property δ ( a b ) = δ ( b ) / | a | d , and therefore Eq. (34) reduces to W ( r , p ; t ) = (cid:90) (cid:90) e imω (cid:126) sin ωt ( u cos( ωt ) − r ) v ρ (cid:16) u + v , u − v (cid:17) δ (cid:18) r cos ωt − u − sin ωtmω p (cid:19) d u d v . After u -integration, we get W ( r , p ; t ) = (cid:90) e − i (cid:126) [ mω r sin ωt + p cos ωt ] . v ρ (cid:18) r cos ωt − sin ωtmω p + v , r cos ωt − sin ωtmω p − v (cid:19) d v . Notice that according to Eq. (7), the above expression represents the initial Wigner function at a multidimensionalphase space point (cid:16) r cos ωt − p sin ωtmω , mω r sin ωt + p cos ωt (cid:17) , that is W ( r , p ; t ) = W (cid:18) r cos ωt − p sin ωtmω , mω r sin ωt + p cos ωt (cid:19) , (36)which ends the proof and reduces to the result given in Eq. (31) for d = 1. Notice that if ω →
0, this relation reducesto the result given in Eq. (18), as it should.The existence of alternative derivations for a given problem can enrich one’s insight in solving other problems relatedto it. In this respect, it is interesting to see how our generalized phase space result translates in a real space. In otherwords, what is the relation between the one-body density matrix ρ ( r , r ; t ) at time t and its initial value ρ ( r , r ) at t = 0 in d dimensions. To obtain this connection, we shall directly apply the inverse Wigner transformation to bothsides of Eq. (8) with (36), we obtain the one-body density matrix in real space at time tρ ( r + s / , r − s / t ) = (cid:90) d p (2 π (cid:126) ) d W (cid:18) r cos ωt − p sin ωtmω , mω r sin ωt + p cos ωt (cid:19) e i p · s / (cid:126) . (37)Although the analytical expression of W in d = 1 is known to be given in terms of Laguerre polynomials L k forarbitrary N [19] as W ( x, p ) = 2 exp (cid:40) − (cid:126) ω (cid:18) p m + 12 mω x (cid:19) N − (cid:88) k =0 ( − k L k (cid:18) (cid:126) ω (cid:18) p m + 12 mω x (cid:19)(cid:19)(cid:41) , and its expression in d dimensions can be shown to be given in terms of generalized Laguerre polynomials L ( d − k [25], W ( r , p ) = 2 d exp (cid:40) − (cid:126) ω (cid:18) p m + 12 mω r (cid:19) N − (cid:88) k =0 ( − k L ( d − k (cid:18) (cid:126) ω (cid:18) p m + 12 mω r (cid:19)(cid:19)(cid:41) , (38)here its use to calculate the integral in Eq. (37) is not straightforward.Alternatively, we propose to use expression of the so-called Bloch propagator in Wigner phase space which is relatedto the density operator through an appropriate Laplace transform. This relation is obtained as follows. For a systemof N noninteracting fermions moving in a potential V ( r ), the one-body density matrix is ρ ( r , r ) = (cid:88) k φ k ( r ) φ ∗ k ( r ) θ ( µ − ε k ) , where the sum is over occupied single-particle states up to the Fermi energy µ , and the φ k ’s and ε k ’s are respectivelythe normalized single particle wavefunctions and their corresponding energies, that is H φ k = ε k φ k . The Heavisideunit-step function is denoted by θ ( x ) and by using the inverse Laplace transform identity [34] θ ( µ − ε k ) = (cid:90) c + i ∞ c − i ∞ dz πi e z ( µ − ε k ) z , we can write the density matrix ρ ( r , r ) as an inverse Laplace transformation (see [18] and references quoted therein),so that ρ ( r , r ) = c + i ∞ (cid:90) c − i ∞ dz πi e zµ C ( r , r ; z ) z , (39)with C ( r , r ; z ) = (cid:80) k φ k ( r ) φ ∗ k ( r ) e − zε k the matrix elements of the Bloch operator e − zH . Here the parameter z isconsidered as a mathematical variable which in general is taken to be complex and c is a positive constant. We canimmediately obtain the desired relation by writing the Wigner phase space version of Eq. (39) so that W ( r , p ) = c + i ∞ (cid:90) c − i ∞ dz πi e zµ (cid:101) C ( r , p ; z ) z , (40)0where (cid:101) C ( r , p ; z ) denotes the Wigner transform of C ( r , r ; z ). For the case of an isotropic harmonic potential in d dimensions, the phase space function (cid:101) C ( r , p ; z ) has a simple explicit expression [17, 19], (cid:101) C ( r , p ; z ) = 1 (cid:2) cosh (cid:0) z (cid:126) ω (cid:1)(cid:3) d e − (cid:126) ω ( tanh z (cid:126) ω ) (cid:16) p m + mω r (cid:17) . (41)To end these preparations, we give the expression of the inverse Wigner transform of Eq. (41), in terms of centre ofmass and relative coordinates (see [35] and references quoted therein), C ( r + s / , r − s /
2) = (cid:16) mω π (cid:126) sinh z (cid:126) ω (cid:17) d e − mω (cid:126) (cid:104) r tanh ( z (cid:126) ω ) + s coth ( z (cid:126) ω ) (cid:105) , (42)an expression which will be used shortly.Now we are in position to proceed with the integral in Eq. (37). Let us substitute Eq. (40) into (37), we obtain ρ ( r + s , r − s ; t ) = c + i ∞ (cid:90) c − i ∞ dz πi e zµ z (cid:90) d p (2 π (cid:126) ) d e i p · s (cid:126) (cid:101) C (cid:18) r cos ωt − p sin ωtmω , mω r sin ωt + p cos ωt ; z (cid:19) , (43)and using the expression of (cid:101) C in Eq. (41), we get ρ ( r + s / , r − s / t ) = c + i ∞ (cid:90) c − i ∞ dz πi e zµ z (cid:2) cosh (cid:0) z (cid:126) ω (cid:1)(cid:3) d (cid:90) d p (2 π (cid:126) ) d e i p · s (cid:126) × (cid:34) e − (cid:126) ω ( tanh z (cid:126) ω ) (cid:18) ( mω r sin ωt + p cos ωt )22 m + mω ( r cos ωt − p sin ωtmω ) (cid:19) (cid:35) . (44)The above integral on p is carried out in the Appendix B, and therefore we have ρ ( r + s , r − s ; t ) = e i m (cid:126) ˙ b ( t ) b ( t ) r · s ( b ( t )) d c + i ∞ (cid:90) c − i ∞ dz πi e zµ z (cid:18) mω π sinh z (cid:126) ω (cid:19) d e − mω (cid:126) (cid:104) r b tanh ( z (cid:126) ω ) + s b coth ( z (cid:126) ω ) (cid:105) , (45)where b ( t ) is a time dependent scaling factor given by b ( t ) = (cid:115) (cid:18) ω ω − (cid:19) sin ωt, (46)which is the solution of ˙ b ( t ) = (cid:16) ω ω − (cid:17) ( ω sin ωt cos ωt ) /b ( t ). Equation (42) allows us to write Eq. (45) as ρ ( r + s / , r − s / t ) = e i m (cid:126) ˙ b ( t ) b ( t ) r · s ( b ( t )) d c + i ∞ (cid:90) c − i ∞ dz πi e zµ z C (cid:18) r b ( t ) + s b ( t ) , r b ( t ) − s b ( t ) (cid:19) . (47)We observe that the above inverse Laplace transform is nothing but the initial one-body density matrix at rescaledpositions, with b ( t ) as the scaling factor. We then arrive to ρ ( r + s / , r − s / t ) = e i m (cid:126) ˙ b ( t ) b ( t ) r · s ( b ( t )) d ρ (cid:18) r b ( t ) + s b ( t ) , r b ( t ) − s b ( t ) (cid:19) , (48)returning to the original coordinates, r = r + s / r = r − s /
2, and since r · s = ( r − r ) /
2, the above scaling lawbecomes ρ ( r , r ; t ) = e i m (cid:126) ˙ b ( t ) b ( t ) ( r − r ) ( b ( t )) d ρ (cid:18) r b ( t ) , r b ( t ) (cid:19) , (49)a result valid for arbitrary dimensions. For d = 1, one recovers the result obtained by the rescaling method [27].1 A. Ballistic versus non-ballistic expansions in phase space
It is important to note that the above scaling law obtained for trap to trap quench is not restricted just tononinteracting particles. In fact, the scaling law in Eq. (49) was proven to hold also for some physical systemsof interacting particles in which the interactions are acting before and after the quench of the harmonic potential(non ballistic expansion). This is the case for the two following situations: (i) the Tonks-Girardeau gas, whichconsists in a gas of identical bosons interacting through very strong repulsive zero-range interactions, confined byharmonic trap in one dimension ( d = 1) [36, 37], and (ii) for harmonically trapped interacting fermions in threedimensions ( d = 3) at unitarity [38]. Notice that long time ago Pitaevskii and Rosch [39] introduced a scaling ansatzfor a two dimensional bosonic system of particles interacting with contact or inverse square interaction. Later on,scaling approach to quantum non-equilibrium dynamics of interacting systems subject to external linear and parabolicpotentials has been examined in [40], where many-body scaling solutions to more general types of interaction andarbitrary dimensionality where obtained.It may be of interest to see how the above scaling law is expressed in Wigner phase space. For obtaining this result,we substitute Eq. (48) into (7) and we obtain W ( r , p ; t ) = 1( b ( t )) d (cid:90) d s ρ (cid:18) r b ( t ) + s b ( t ) , r b ( t ) − s b ( t ) (cid:19) e − i (cid:126) ( b ( t ) p − m ˙ b ( t ) r ) · s b ( t ) = (cid:90) d u ρ (cid:18) r b ( t ) + u , r b ( t ) − u (cid:19) e − i (cid:126) ( b ( t ) p − m ˙ b ( t ) r ) · u , (50)and according to Eq. (7), the right-hand side represent the initial Wigner function at phase space point (cid:16) r b ( t ) , b ( t ) p − m ˙ b ( t ) r (cid:17) , so that W ( r , p ; t ) = W (cid:18) r b ( t ) , b ( t ) p − m ˙ b ( t ) r (cid:19) , (51)where we recall that b ( t ) is given by (46).As stated for Eq. (49), it follows that the relation in Eq. (51) is not only valid for noninteracting particle systemsbut also for the above two systems pertaining to interacting particles. The interested reader may ask on the differencebetween Wigner functions given respectively by Eq. (36) and (51). The time dependent Wigner function in Eq. (36)describes a ballistic expansion (interactions are suppressed) following the quench of the harmonic potential whilethe Wigner function in Eq. (51) concerns ballistic or a nonballistic expansion after the quench of the harmonic trappertaining to the nature (interacting or noninteracting) of the initially confined system before the quench. We canshow that for an initially harmonically confined noninteracting particle system subjected to a quench of the potential,the two forms Eq. (36) and (51) are identical, that is W (cid:18) r cos ωt − p sin ωtmω , mω r sin ωt + p cos ωt (cid:19) = W (cid:18) r b ( t ) , b p − m ˙ b ( t ) r (cid:19) . (52)In fact for a system of N noninteracting fermions confined in a d dimensional isotropic harmonic potential withfrequency ω , the Wigner function was given in (38), where we can observe that W ( r , p ) depends on the phase spacevariables ( r , p ) by means of the classical Hamiltonian H cl ( r , p ) = (cid:16) p m + mω r (cid:17) , and is nothing but the Wignertransform of the quantum one particle Hamiltonian. Hence, to prove Eq. (52) one has just to check the equality( mω r sin ωt + p cos ωt ) m + mω (cid:18) r cos ωt − p sin ωtmω (cid:19) = (cid:16) b ( t ) p − m ˙ b ( t ) r (cid:17) m + mω (cid:18) r b ( t ) (cid:19) . Recalling (46) it is easy to verify the last equation.We believe that for interacting particles Eq. (52) is no longer true, therefore in that case W (cid:18) r cos ωt − p sin ωtmω , mω r sin ωt + p cos ωt (cid:19) (cid:54) = W (cid:18) r b ( t ) , b ( t ) p − m ˙ b ( t ) r (cid:19) . As stated before, Eq. (36) is valid for both interacting or noninteracting system of particles providing that the gasexpands ballistically (suppression of interactions) after the quench. In this respect we shall exploit this relation to seehow one can access to the initial momentum density of the interacting system using the expansion of the density inreal space.2
B. Recovering the initial momentum density
Recently, an experimental technique to directly image the momentum distribution of a strongly interacting two-dimensional quantum gas was obtained and characterized [41]. This method is based on the fact that, just afterswitching-off the initially confining trap, and instead of a free expansion, the gas is subjected to an external harmonicpotential V ( r ) = mω r / T = 2 π/ω , the spatial distribution is related to the momentum density ofthe initially confined quantum gas. In the following, we provide a generalization of this relation in arbitrary spatialdimensions by exploiting the relation in Eq. (36) for the Wigner function.Let us denote by T the period corresponding to frequency ω of the harmonic trap and considering the specific time t = T / π/ ω after the quench, Eq. (36) reduces to W (cid:18) r , p ; T (cid:19) = W (cid:16) − p mω , mω r (cid:17) . Using Eq. (9), the spatial density at this time is ρ (cid:18) r ; T (cid:19) = (cid:90) d p (2 π (cid:126) ) d W (cid:16) − p mω , mω r (cid:17) . Making the change of variable u = − p / ( mω ), the above integral in d dimensions becomes ρ (cid:18) r ; T (cid:19) = ( mω ) d (cid:90) d u (2 π (cid:126) ) d W ( u , mω r ) , (53)and according to Eq. (10), the right-hand side of Eq. (53) essentially represents the initial momentum density, that is ρ (cid:18) r ; T (cid:19) = ( mω ) d n ( p = mω r ; 0) , (54)a relation which clearly exhibits the mapping between the initial momentum density for a given value of the momentum p = mω r and the spatial density at position r at time t = T / ω . IV. CONCLUDING REMARKS
In this paper we have studied the non-equilibrium dynamics in phase space generated by a sudden change of theHamiltonian in a quantum system, through the analysis of the Wigner function. For the case of two attractive particleswe calculated the corresponding time dependent Wigner function following a swich-off of the interaction.For possible experimental implementation in ultra-cold quantum gases field, an interaction of the form − aδ ( x ) + bδ (cid:48) ( x ) , a > , b ∈ R was considered and we have calculated the two-particle Wigner function. At large values ofmomentum p , we have found that the associated momentum distribution scales as 1 /p , while in the absence of δ (cid:48) interaction, this density scales, in this case of pure zero-range delta interaction, as the well known law 1 /p .We have generalized to arbitrary dimensions d a derivation, by using alternative method, of a relationship shownrecently in one dimension between the Wigner function at time t and its initial value following a sudden change ofthe harmonic trap of noninteracting particles. We have exploited our generalized relation, through the use of inverseWigner transformation to obtain in d dimensions the scaling law satisfied by the one-body density matrix in realspace. Using the generalized relation in Wigner phase space for the considered quench, we have shown that the initialmomentum density of a system of particles (interacting or noninteracting) is exactly mapped for a given value of themomentum p = mω r to the spatial density at position r at time a quarter of time period, t = T /
4, after the ballisticexpansion in the harmonic trap with frequency ω . It should be noted that our method can be easily adapted to dealwith physical situations corresponding to quasi-1 d or quasi-2 d configurations.An interesting extension of this dynamically situation would be to study the problem of a Lieb-Liniger gas at finiterepulsion strength. Work in this direction is in progress. Acknowledgements
This work was partially supported by the Spanish MINECO (MTM2014-57129-C2-1-P), Junta de Castilla y Le´onand FEDER projects (BU229P18, VA057U16, and VA137G18). K.B. thanks the Direction G´en´erale de la Recherche3Scientifique et du D´eveloppement Technologique (DGRSDT-Algeria) for financial support. The authors acknowledgethe anonymous referee for helpful suggestions.
Appendix A
In this Appendix we shall prove the relationship given in Eq. (6). Using Dirac notations, let us rewrite Eq. (5) as ρ ( r , r (cid:48) ; t ) = N (cid:90) d r · · · d r N (cid:104) r , r , . . . , r N | Φ( t ) (cid:105) (cid:104) Φ( t ) | r (cid:48) , r , . . . , r N (cid:105) = N (cid:90) d r · · · d r N (cid:104) r , r , . . . , r N | e − i (cid:126) Ht | Φ (cid:105) (cid:104) Φ | e + i (cid:126) Ht | r (cid:48) , r , . . . , r N (cid:105) , (55)where we have used Eq. (4) to obtain the second form. Applying the closure relationship, we can write ρ ( r , r (cid:48) ; t ) = N (cid:90) d r · · · d r N (cid:90) d r (cid:48) · · · d r (cid:48) N (cid:90) d r (cid:48)(cid:48) · · · d r (cid:48)(cid:48) N Φ ( r (cid:48) , . . . , r (cid:48) N )Φ ∗ ( r (cid:48)(cid:48) , . . . , r (cid:48)(cid:48) N ) (56) × (cid:104) r , r , . . . , r N | e − i (cid:126) Ht | r (cid:48) , r (cid:48) , . . . , r (cid:48) N (cid:105) (cid:104) r (cid:48)(cid:48) , r (cid:48)(cid:48) , . . . , r (cid:48)(cid:48) N | e + i (cid:126) Ht | r (cid:48) , r , r , . . . , r N (cid:105) . Since the post-quench Hamiltonian H describes a system of N noninteracting particles, we can write the followingfactorization for the time evolution operator, so that e − i (cid:126) Ht = e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) · · · e − it (cid:126) (cid:18) p N m + V ( r N ) (cid:19) , (57)from which we immediately deduce (cid:104) r , r , · · · , r N | e − it (cid:126) H | r (cid:48) , r (cid:48) , . . . , r (cid:48) N (cid:105) = (cid:104) r | e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) · · · (cid:104) r N | e − it (cid:126) (cid:18) p N m + V ( r N ) (cid:19) | r (cid:48) N (cid:105) , and similarly (cid:104) r (cid:48)(cid:48) , r (cid:48)(cid:48) , . . . , r (cid:48)(cid:48) N | e it (cid:126) H | r (cid:48) , r , r , . . . , r N (cid:105) = (cid:104) r (cid:48)(cid:48) | e it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) · · · (cid:104) r (cid:48)(cid:48) N | e it (cid:126) (cid:18) p N m + V ( r N ) (cid:19) | r N (cid:105) . Substituting these last two relations into Eq. (56) we obtain ρ ( r , r (cid:48) ; t ) = N (cid:90) d r · · · d r N (cid:90) d r (cid:48) · · · d r (cid:48) N (cid:90) d r (cid:48)(cid:48) · · · d r (cid:48)(cid:48) N Φ ( r (cid:48) , . . . , r (cid:48) N )Φ ∗ ( r (cid:48)(cid:48) , , . . . , r (cid:48)(cid:48) N ) × (cid:104) r | e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) (cid:104) r | e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) · · · (cid:104) r N | e − it (cid:126) (cid:18) p N m + V ( r N ) (cid:19) | r (cid:48) N (cid:105)× (cid:104) r (cid:48)(cid:48) | e it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) (cid:104) r (cid:48)(cid:48) | e it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:105) · · · (cid:104) r (cid:48)(cid:48) N | e it (cid:126) (cid:18) p N m + V ( r N ) (cid:19) | r N (cid:105) , which can be rewritten as follows ρ ( r , r (cid:48) ; t ) = N (cid:90) d r (cid:48) · · · d r (cid:48) N (cid:90) d r (cid:48)(cid:48) · · · d r (cid:48)(cid:48) N Φ ( r (cid:48) , . . . , r (cid:48) N )Φ ∗ ( r (cid:48)(cid:48) , . . . , r (cid:48)(cid:48) N ) × (cid:104) r | e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) (cid:104) r (cid:48)(cid:48) | e it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) (cid:90) d r (cid:104) r (cid:48)(cid:48) | e it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:105) (cid:104) r | e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) × · · · (cid:90) d r N (cid:104) r (cid:48)(cid:48) N | e it (cid:126) (cid:18) p N m + V ( r N ) (cid:19) | r N (cid:105) (cid:104) r N | e − it (cid:126) (cid:18) p N m + V ( r N ) (cid:19) | r (cid:48) N (cid:105) . (58)Using the closure relations over the kets | r (cid:105) , | r (cid:105) . . . . | r N (cid:105) and then carrying out the integrals over the variables r (cid:48)(cid:48) , . . . , r (cid:48)(cid:48) N , the above expression reduces to ρ ( r , r (cid:48) ; t ) = N (cid:90) (cid:90) d r (cid:48) d r (cid:48)(cid:48) (cid:104) r | e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) (cid:104) r (cid:48)(cid:48) | e it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105)× (cid:90) d r (cid:48) · · · d r (cid:48) N Φ ( r (cid:48) , r (cid:48) , . . . , r (cid:48) N )Φ ∗ ( r (cid:48)(cid:48) , r (cid:48) , . . . , r (cid:48) N ) . (59)4Taking into account the definition of the initial reduced one-body density matrix in Eq. (2), we arrive to ρ ( r , r (cid:48) ; t ) = N (cid:90) (cid:90) d r (cid:48) d r (cid:48)(cid:48) (cid:104) r | e − it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) (cid:104) r (cid:48)(cid:48) | e it (cid:126) (cid:18) p m + V ( r ) (cid:19) | r (cid:48) (cid:105) ρ ( r (cid:48) , r (cid:48)(cid:48) ) (60)which we rewrite as ρ ( r , r (cid:48) ; t ) = (cid:90) (cid:90) d r (cid:48) d r (cid:48)(cid:48) U ( r , r (cid:48) ; t ) U ∗ ( r (cid:48) , r (cid:48)(cid:48) ; t ) ρ ( r (cid:48) , r (cid:48)(cid:48) ) , (61)where U ( r , r (cid:48) ; t ) = (cid:104) r | e − it (cid:126) ( p m + V ( r )) | r (cid:48) (cid:105) is the single particle propagator associated to the many-body noninteractingpost-quench Hamiltonian. Now changing the names of the variables in Eq. (61), so that, r (cid:48) = ξ , r (cid:48)(cid:48) = ξ , the desiredresult in Eq. (6) is recovered. Appendix B. Derivation of Eq. (45)
To calculate the integral on p in Eq. (44), we need some intermediate results. First of all, it is easy to check that( mω r sin ωt + p cos ωt ) m + 12 mω (cid:18) r cos ωt − p sin ωtmω (cid:19) = (62)= 12 mω (cid:20) ω ω − b ( t ) (cid:21) r + 12 m b ( t ) p − ω (cid:18) ω ω − (cid:19) (sin ωt cos ωt ) r · p where b ( t ) = (cid:114) (cid:16) ω ω − (cid:17) sin ωt is a time dependent scaling factor. After insertion of this result in Eq. (44), itbecomes ρ ( r + s / , r − s / t ) = c + i ∞ (cid:90) c − i ∞ dz πi e zµ z (cid:2) cosh (cid:0) z (cid:126) ω (cid:1)(cid:3) d e − mω (cid:126) ω ( tanh z (cid:126) ω ) (cid:20) ω ω − b ( t ) (cid:21) r (63) × (cid:90) d p (2 π (cid:126) ) d e i p · s (cid:126) e − (cid:126) ω ( tanh z (cid:126) ω ) (cid:18) m b ( t ) p − ω (cid:18) ω ω − (cid:19) (sin ωt cos ωt ) r · p (cid:19) . Using the d -dimensional identity (cid:90) d u e − α u e − a · u = (cid:16) πα (cid:17) d e a α , (64)the integration on p in Eq. (63) can then be performed, yielding to the expression of the density matrix in the form ρ ( r + s , r − s ; t ) = 1(2 π (cid:126) b ( t )) d c + i ∞ (cid:90) c − i ∞ dz πi (cid:34) π (cid:126) mω tanh z (cid:126) ω (cid:35) d e zµ e − mω (cid:126) ω ( tanh z (cid:126) ω ) (cid:20) ω ω − ( b ( t )) ( t ) (cid:21) r z (cid:2) cosh (cid:0) z (cid:126) ω (cid:1)(cid:3) d × e mω (cid:126) ω b ( t ))2 ( tanh z (cid:126) ω ) (cid:18) ω ω − (cid:19) (sin ωt cos ωt ) r − mω (cid:126) ( b ( t ))2 ( coth z (cid:126) ω ) s + i mω (cid:126) b ( t ))2 (cid:18) ω ω − (cid:19) (sin ωt cos ωt ) r · s , which after collecting similar terms and simplification turns into Eq. (45). [1] S. Giorgini, L.P. Pitaevski, S. Stringari, Rev. Mod. Phys. , 1215 (2008).[2] I. Bloch, J. Dalibard, W. Zwerger, Rev. Mod. Phys. , 885 (2008).[3] L. D’Alessio, Y. Kafri, A. Polkovnikov, M. Rigol, Adv. Phys. , 239 (2016).[4] F. Borgonovi, F.M. Izrailev, L.F. Santos, V.G. Zelevinsky, Phys. Rep. , 1 (2016).[5] M. Rigol, V. Dunjko, V. Yurovsky, M. Olshanii, Phys. Rev. Lett. , 050405 (2007).[6] M. Rigol, V. Dunjko, M. Olshanii, Nature , 854 (2008). [7] J.S. Caux and R. M. Konik, Phys. Rev. Lett. , 175301 (2012).[8] T. Mori and N. Shiraishi, Phys. Rev. E. , 022153 (2017).[9] M. Eckstein, M. Kollar, P. Werner, Phys. Rev. Lett. , 056403 (2009).[10] J. Cardy, Phys. Rev. Lett. , 220401 (2014).[11] M. Collura, S. Sotiriadis, P. Calabrese, Phys. Rev. Lett. , 245301 (2013).[12] R.M. Dreizler and E.K.U. Gross, Density Functional Theory (Springer, Berlin, 1990).[13] E.P. Wigner, Phys. Rev. , 749 (1932).[14] H.J. Groenewold, Physica , 405 (1946).[15] J.E. Moyal, Proc. Cambridge Philos. Soc. , 99 (1949).[16] M. Gadella, Fortschritte der Physik , 229 (1995).[17] A.M. Ozorio de Almeida and J.H. Hannay, Ann. Phys. , 115 (1982).[18] M. Brack, R.K. Bhaduri, Semiclassical Physics , Frontiers in Physics, vol. 96 (Westview, Boulder, 2003)[19] M. Hillery, R.F. O’Connel, M.O. Scully, E.P. Wigner, Phys. Rep. , 121 (1984).[20] D.F. Walls and G.J. Milburn,
Quantum Optics (Springer, Berlin, 2008).[21] T. Douce et al. , Scientific Reports , 3530 (2013).[22] M.V. Berry, Phil. Trans. R. Soc. London A , 237 (1977).[23] E. Bettelheim and P.B. Wiegmann, Phys. Rev. B. , 085102 (2011).[24] D.S. Dean, P. Le Doussal, S.N. Majumdar, G. Schehr, Phys. Rev. A , 063614 (2018).[25] B.P. van Zyl, J. Phys. A: Math. Theor. , 315302 (2012).[26] K. Bencheikh and L. M. Nieto, Phys. Rev. A , 053614 (2008).[27] D. S. Dean et al, EPL Quantum Mechanics and Path Integrals (McGraw-Hill, New York, 1964).[29] M. Olshanii and V. Dunjko, Phys. Rev. Lett. , 090401 (2003).[30] M. Gadella, J. Negro and L.M. Nieto, Phys. Lett. A , 1310 (2009).[31] M. Gadella, M.L. Glasser and L.M. Nieto, Int. J. Theor. Phys. , 015204 (2016).[33] J.M. Mu˜noz-Casta˜neda, L.M. Nieto, and C. Romaniega, Ann. Phys. , 246 (2018).[34] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1972)[35] K. Bencheikh and L. M. Nieto, J. Phys. A: Math. Theor. (2010) 375002.[36] A. Minguzzi and D.M. Gangardt, Phys. Rev. Lett. , 240404 (2005).[37] P. Ruggiero, Y. Brun, and J. Dubail, SciPost Phys. , 051 (2019).[38] Y. Castin, C. R. Physique , 407 (2004).[39] P. Pitaevskii and A. Rosch, Phys. Rev. A 55, R 853 (1997)[40] V. Gritsev, P. Barmettler and E. Demler, New Journal of Physics 12 (2010) 113005[41] P.A. Murthy et al. , Phys. Rev. A90