Exact results for nonequilibrium dynamics in Wigner phase space

Time evolution of the Wigner distribution function of an initially interacting one-dimensional quantum gas following an interaction quench is examined. Considering the scenario of a quench at $t=0$ from a nonzero to zero (free particles) value of the interaction strength, we derive a simple relationship between the dynamical Wigner distribution function and its initial value. As an application, a two-particle system initially interacting through two different zero-range interactions of Dirac delta type is presented. For a system of particles (interacting or noninteracting) that is suddenly let to move ballistically in a harmonic trap in $d$ dimensions, we derive a relationship between the resulting time dependent Wigner function and its initial value. Our method is based on the use of the time evolution of the reduced one-body density matrix. By means of an inverse Wigner transform we obtain, for the case of an initially harmonically trapped noninteracting particle system in $d$ dimensions, the scaling law satisfied by the time dependent density at time $t$ after a sudden change of the frequency of trapping potential. Ballistic versus nonballistic expansions are analyzed in phase space for physical systems of particles with dynamics are governed by scaling laws.

Time evolution of the Wigner distribution function of an initially interacting one-dimensional quantum gas following an interaction quench is examined. Considering the scenario of a quench at t = 0 from a nonzero to zero (free particles) value of the interaction strength, we derive a simple relationship between the dynamical Wigner distribution function and its initial value. As an application, a two-particle system initially interacting through two different zero-range interactions of Dirac delta type is presented.
For a system of particles (interacting or noninteracting) that is suddenly let to move ballistically in a harmonic trap in d dimensions, we derive a relationship between the resulting time dependent Wigner function and its initial value. Our method is based on the use of the time evolution of the reduced one-body density matrix. By means of an inverse Wigner transform we obtain, for the case of an initially harmonically trapped noninteracting particle system in d dimensions, the scaling law satisfied by the time dependent density at time t after a sudden change of the frequency of trapping potential. Ballistic versus nonballistic expansions are analyzed in phase space for physical systems of particles with dynamics are governed by scaling laws.

I. INTRODUCTION
The description of time evolution properties of quantum many body systems is presently a fundamental topic of research. The amazing development of trapping and cooling techniques have led to experimental realization of quantum many-body systems consisting of ultracold atomic gases where atoms are confined in traps [1,2]. These artificial many-body systems can be produced and loaded in various geometric traps. Such experimental progress allowed a full control of the external parameters in the Hamiltonian governing the quantum system dynamics. An interesting issue in the field of ultracold quantum gases is to study the time evolution of a non-equilibrium situation generated through a quantum quench, which consists of a sudden change of the Hamiltonian parameters (for example a change of the harmonic trap frequency or a change in the interaction strength between the atoms of the gas through Feshbach resonance). A quantum quench is the easiest way to drive a system to non-equilibrium: the system is supposed to be in its Hamiltonian ground state until time t = 0, when the sudden change of a coupling leads to a new Hamiltonian according to which the system evolves for t > 0. On the theoretical side, significant advances have been carried 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 new ensemble in Statistical Mechanics called generalized Gibbs ensemble, which is a more general concept than the usual grand-canonical ensemble and turns out to be a powerful tool in the prediction of relaxation processes for certain integrable one-dimensional systems [5][6][7][8][9][10][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 underlying Hamiltonians before and after the interaction quench. For t < 0 the gas is in equilibrium and its many-body Hamiltonian is where p i = −i ∇ i is the momentum operator of the particle i, V 0 is an external confining potential, and v(r i , r j ) is a two-body atom-atom interaction. We assume that the system is initially in a quantum many-body state |Φ 0 and its associated initial reduced one-body-density matrix ρ 0 (r, r ) is defined as (see for instance Ref. [12]) At t = 0, the interactions are turned-off and the many-body Hamiltonian is given by 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 It should be noted that |Φ 0 is not an eigenstate of the post-quench Hamiltonian H. The latter Hamiltonian describes a system of independent particles, while |Φ 0 describes the state of an initially interacting particles. Before proceeding further, we shall first derive a relationship between the initial ρ 0 (r, r ), and the time dependent ρ(r, r ; t) reduced one-body-density matrices. Using the definition we show in Appendix A the following relation where U (r, ξ 1 ; t) = r| e − i ( p 2 2m +V )t |ξ 1 is the single particle propagator associated to the post quench Hamiltonian in Eq. (3). The above relation describes the time evolution of the one-body density matrix for the considered specific quench 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 this quench scenario we relate the so-called dynamical Wigner distribution function to its initial value just before the quench. Our interest in the Wigner distribution function [13] is motivated by the fact that it provides a useful tool to study various properties of many-body systems. Besides, it is well known that it allows a reformulation of quantum mechanics in terms of classical concepts [14,15], and it is also used to generate semi-classical approximations [16,17]. Although other distribution functions exist, the Wigner distribution function has the virtue of its mathematical simplicity. Nevertheless, it may take negative values as a manifestation of its quantum nature, and therefore it does not represent a true probability but a quasi-probability distribution [16,18]. Wigner distribution functions have been used in various contexts, as in cold atomic gases [1,2], quantum optics [19], quantum information [20], quantum chaos [21], and in the study of non-equilibrium dynamics generated by the perturbation of a Fermi gas system [22]. Also, the Wigner distribution function of the noninteracting limit of a Fermi gas system at zero and nonzero temperatures has been the subject of recent studies [23][24][25]. The Wigner function is defined as the Fourier transform of the one-body density matrix ρ(r 1 , r 2 ; t) ≡ ρ(r + s/2, r − s/2; t) on the relative coordinate s = r 1 − r 2 , where r = (r 1 + r 2 )/2 is the centre of mass coordinate. To be precise, we have in d dimensions W (r, p; t) = ds ρ(r + s/2, r − s/2; t) e −ip.s/ , which is a function of phase space variables (r, p) at time t. The inverse transformation reads ρ(r + s/2, r − s/2; t) = dp (2π ) d W (r, p; t)e +ip.s/ .
The structure of the paper is as follows. In Section II we consider the case of an initially untrapped one dimensional interacting system subjected to a sudden swich-off of interactions. The resulting ballistic dynamics is examined in phase space. We derive a relationship between the dynamical Wigner phase space density at time t > 0 and its initial value before the quench. As an application, we analyze the case of a two-particle system initially interacting through two different zero-range interactions of Dirac delta type. Section III is devoted to a study of nonequilibrium dynamics through the Wigner function of harmonically trapped noninteracting particles. Very recently, this system, in one spatial dimension (d = 1), has been the subject of an interesting study [26], where it is suddenly subjected to a modification of the trapping frequency and a relationship between the resulting time dependent Wigner function and its initial value is obtained. In this section we propose a generalization to arbitrary spatial dimensions (d ≥ 1) of this relation by using an alternative fully quantum mechanical method which is based on the use of time evolution of the one-body density matrix. As a bonus, by using the inverse Wigner transformation to our generalized relationship in Wigner 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 obtained in [26], as it should.
In order to observe how the dynamical Wigner function following a quench is affected by interactions, we consider physical systems of particles whose dynamics are governed by scaling laws suffering a non-ballistic expansion, and we compare the resulting Wigner function with the one obtained for a ballistic expansion. Some final conclusions put an end to the paper in Section IV.

II. TIME EVOLUTION OF WIGNER FUNCTION OF AN INITIALLY UNTRAPPED INTERACTING SYSTEM FOLLOWING A SUDDEN SWICH-OFF OF INTERACTIONS.
In this section let us consider the situation where both external potentials V 0 and V in Eqs. (1) and (3) vanish, and the study is restricted to one-dimensional interacting particles system. The time evolution of the resulting reduced one-body density matrix after a finite time of free expansion t > 0 is then given by the one dimensional version of Eq. (6), that is Here The matrix element U (x, ξ; t) is the free particle Feynman propagator in the configuration space, given by [27] U If we substitute this expression into Eq. (11) we find which can be rewritten as Now it is more convenient to introduce the center-of-mass and relative coordinates, respectively defined by z = (ξ 1 + ξ 2 )/2 and z = ξ 1 − ξ 2 . Then we can write 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 Carrying out the integration on the variable s, we obtain which after integration over z reduces to But according to Eq. (7), the above integral is nothing but the initial Wigner distribution function at t = 0, and at the phase space point (x − pt/m, p), so that for t > 0 This relation is the Wigner phase space version of Eq. (11). Upon integration over x and making use of the definition in Eq. (10), we get for t > 0 Therefore, we find that under ballistic expansion (where interatomic collisions during the expansion are not present) the dynamical momentum density for t > 0 is equal to its initial value n(p, 0). Some information can be obtained from this result. For example, if one considers an initial quantum gas with short range interaction, whose momentum distribution n(p, 0) exhibits a 1/p 4 tail at large momentum (see for instance Ref. [28]), equation (19) says that this long 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 0 = 0) provided that one realizes at t = 0 a simultaneously sudden double quench, where the interactions are turned-off with the release 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 Dirac delta potential with highly asymmetric mass imbalance: the particle with mass M is so heavy that the center-of-mass motion can be ignored. The problem reduces then to a one-body system and the m mass light particle Hamiltonian is given for t ≤ 0 by Here a > 0 is the strength of the interaction. The single normalized bound state is φ(x) = √ α e −α|x| with energy E a = −ma 2 /(2 2 ), and α = ma/ 2 . The corresponding Wigner distribution function for t ≤ 0 can be calculated analytically using (7) and is found to be 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 Wigner function for t > 0: A plot of this function W (x, p; t) is given for a particular value of t > 0 on the right hand side of Figure 1. A clear distortion can be observed as time grows: for t = 0 there are two symmetry axes, x = 0 and p = 0; for bigger values of t the perfect symmetry is lost, although the axes p = 0 and p = mx/t (this one clockwise rotating with time) play an important role; squeezing is more and more pronounced as t → ∞, the axis p = mx/t approaching the axis p = 0. Finally, using (19) it is easy to show analytically that, the initial momentum density is It has been shown in [28] that for a system of zero-range δ−interacting one-dimensional atoms with arbitrary strength, the high-p asymptotic behavior of the momentum distribution for both free and harmonically trapped atoms, exhibits a universal 1/p 4 dependence. As can be seen in Eq. (23) and at large values of the momentum p, we recovered this 1/p 4 dependence of the momentum distribution.

B. Case of two particles interacting through a δ-δ interaction
We are going to consider now an extension of the previous study of two interacting particles that takes into account the presence of an extra point-like interaction term in the potential, proportional to δ . This type of point or zerorange potentials are a subject of recent study in differents contexts [29][30][31][32]. The Hamiltonian of the light particle with mass m is now given for t < 0 by The associated Schrödinger equation has been carefully analyzed in [29], where it was proven that the above Hamiltonian supports only one bound state of energy The normalized wave function is where and θ(x) stands for the Heaviside unit step function. In this case the Wigner distribution function for t ≤ 0 can be also determined analytically from (7) and it turns out to be the following expression 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. Some plots of this function are given in Figure 2 for At time t = 0 the interaction is turned off (a = b = 0) and then, from (18) we get the following Wigner function for t > 0: Plots of this function W (x, p; t) are given in Figure 3 for t = 1 and values b = −0.9 (left) and b = −1.2 (right). A clear distortion can be observed comparing with Figure 2: the axis p = 0 is preserved, but the axis x = 0 rotates around the origin and becomes x = p (this one changing with time); squeezing is stronger as t → ∞, the axis x = 0 approaching the axis p = 0. From (19) and (28) the initial momentum density can be determined analytically: which coincides with (23) in the limit B = mb/ 2 → 0 and which scales as 1/p 2 for high p if B = 0 (if B = 0 it was already mentioned after (23) that it scales as 1/p 4 ). Remark that in order to have a physical meaning this density n(p, 0) must be a positive function, and therefore we have to impose the following restriction on the coefficient of the δ term in the potential: − √ 2 ≤ B ≤ √ 2, a condition not observed before in the literature. It is interesting to observe that the presence of the additional δ interaction term, changes substantially the high momentum tail of the momentum density. In fact, at large values of the momentum p, the momentum distribution in Eq. (30) exhibits 1/p 2 behavior while in the absence of δ interaction this density scales as 1/p 4 .

III. TIME EVOLUTION OF WIGNER DISTRIBUTION FUNCTION FOLLOWING A SUDDEN CHANGE OF TRAPPING HARMONIC POTENTIAL IN D DIMENSIONS
Very recently Dean et al. [26] computed the time evolution of the Wigner function W (x, p; t) for a one-dimensional system of noninteracting fermions (called also system of independent fermions), initially trapped by a harmonic oscillator potential with frequency ω 0 , subjected to a sudden change of frequency passing from ω 0 to ω. At time t, after this quench it was found (see Eq. (13) of [26]) that where W 0 is the initial Wigner function of the system. To derive this relationship the Liouville equation is used. Here we propose an alternative fully quantum mechanical derivation based on the use of time evolution of one-body density matrix which is valid for an arbitrary spatial dimension d, which obviously reduces to the result in (31) for d = 1 as it should.
Consider a system of noninteracting fermions initially confined in a d−dimensional isotropic harmonic oscillator potential of the form V 0 (r) = mω 2 0 r 2 /2. Suddenly the frequency is changed from ω 0 to ω and the system expands in the potential V (r) = mω 2 r 2 /2. According to (6), the time evolution of the resulting one-body density matrix at t > 0 is then given by Here ρ 0 (ξ 1 , ξ 2 ) is the one-body density matrix at t = 0 and U osc (r 1 , ξ 1 ; t) is the well known propagator associated to the isotropic harmonic oscillator potential with frequency ω, given in d dimensions by [27] U osc (r 1 , ξ 1 ; t) = mω 2πi sin ωt by substitution we obtain Using the centre of mass u = (ξ 1 + ξ 2 )/2 and relative v = (ξ 1 − ξ 2 ) coordinates we can write ρ(r 1 , r 2 ; t) = mω 2π |sin ωt| with r = (r 1 + r 2 )/2 and s = r 1 − r 2 . We can perform the d dimensional s integration, where we have used the property δ(ab) = δ(b)/ |a| d , and therefore Eq. (34) reduces to After u-integration, we get Notice that according to Eq. (7), the above expression represents the initial Wigner function at a multidimensional phase space point r cos ωt − p sin ωt mω , mωr sin ωt + p cos ωt , that is W (r, p; t) = W 0 r cos ωt − p sin ωt mω , mωr sin ωt + p cos ωt , which ends the proof and reduces to the result given in Eq. (31) for d = 1. Notice that if ω → 0, this relation reduces to 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 related to it. In this respect, it is interesting to see how our generalized phase space result translates in a real space. In other words, what is the relation between the one-body density matrix ρ(r 1 , r 2 ; t) at time t and its initial value ρ 0 (r 1 , r 2 ) at t = 0 in d dimensions. To obtain this connection, we shall directly apply the inverse Wigner transformation to both sides of Eq. (8) with (36), we obtain the one-body density matrix in real space at time t ρ(r + s/2, r − s/2; t) = dp (2π ) d W 0 r cos ωt − p sin ωt mω , mωr sin ωt + p cos ωt e ip·s/ .
Although the analytical expression of W 0 in d = 1 is known to be given in terms of Laguerre polynomials L k for arbitrary N [18] as and its expression in d dimensions can be shown to be given in terms of generalized Laguerre polynomials L (d−1) k [24], 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 related to the density operator through an appropriate Laplace transform. This relation is obtained as follows. For a system of N noninteracting fermions moving in a potential V (r), the one-body density matrix is where the sum is over occupied single-particle states up to the Fermi energy µ, and the φ k 's and ε k 's are respectively the normalized single particle wavefunctions and their corresponding energies, that is H 0 φ k = ε k φ k . The Heaviside unit-step function is denoted by θ(x) and by using the inverse Laplace transform identity [33] we can write the density matrix ρ 0 (r 1 , r 2 ) as an inverse Laplace transformation (see [17] and references quoted therein), so that with C 0 (r 1 , r 2 ; z) = k φ k (r 1 )φ * k (r 2 )e −zε k the matrix elements of the Bloch operator e −zH0 . Here the parameter z is considered as a mathematical variable which in general is taken to be complex and c is a positive constant. We can immediately obtain the desired relation by writing the Wigner phase space version of Eq. (39) so that where C 0 (r, p;z) denotes the Wigner transform of C 0 (r 1 , r 2 ; z). For the case of an isotropic harmonic potential in d dimensions, the phase space function C 0 (r, p;z) has a simple explicit expression [16,18], To end these preparations, we give the expression of the inverse Wigner transform of Eq. (41), in terms of centre of mass and relative coordinates (see [34] and references quoted therein), 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 e zµ z dp (2π ) d e i p·s C 0 r cos ωt − p sin ωt mω , mωr sin ωt + p cos ωt; z , The above integral on p is carried out in the Appendix B, and therefore we have where b(t) is a time dependent scaling factor given by which is the solution ofḃ(t) = We observe that the above inverse Laplace transform is nothing but the initial one-body density matrix at rescaled positions, with b(t) as the scaling factor. We then arrive to returning to the original coordinates, r 1 = r + s/2, r 2 = r − s/2, and since r · s = (r 2 1 − r 2 2 )/2, the above scaling law becomes a result valid for arbitrary dimensions. For d = 1, one recovers the result obtained by the rescaling method [26].

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 to noninteracting particles. In fact, the scaling law in Eq. (49) was proven to hold also for some physical systems of 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, which consists in a gas of identical bosons interacting through very strong repulsive zero-range interactions, confined by harmonic trap in one dimension (d = 1) [35], and (ii) for harmonically trapped interacting fermions in three dimensions (d = 3) at unitarity [36].
It is 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 and according to Eq. (7), the right-hand side represent the initial Wigner function at phase space point 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 systems but also for the above two systems pertaining to interacting particles. The interested reader may ask on the difference between 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 while the Wigner function in Eq. (51) concerns ballistic or a nonballistic expansion after the quench of the harmonic trap pertaining to the nature (interacting or noninteracting) of the initially confined system before the quench. We can show 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 , bp−mḃ(t)r .
In fact for a system of N noninteracting fermions confined in a d dimensional isotropic harmonic potential with frequency ω 0 , the Wigner function was given in (38), where we can observe that W 0 (r, p) depends on the phase space variables (r, p) by means of the classical Hamiltonian H cl (r, p) = p 2 2m + 1 2 mω 2 0 r 2 , and is nothing but the Wigner transform 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) 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 0 r cos ωt − p sin ωt mω , mωr sin ωt + p cos ωt = W 0 r b(t) , b(t)p−mḃ(t)r .
As stated before, Eq. (36) is valid for both interacting or noninteracting system of particles providing that the gas expands ballistically (suppression of interactions) after the quench. In this respect we shall exploit this relation to see how one can access to the initial momentum density of the interacting system using the expansion of the density in real space.

B. Recovering the initial momentum density
Recently, an experimental technique to directly image the momentum distribution of a strongly interacting twodimensional quantum gas was obtained and characterized [37]. This method is based on the fact that, just after switching-off the initially confining trap, and instead of a free expansion, the gas is subjected to an external harmonic potential V (r) = mω 2 r 2 /2 where the gas moves ballistically (the interactions are suppressed). It was shown that after a quarter of the oscillator time period T = 2π/ω, the spatial distribution is related to the momentum density of the initially confined quantum gas. In the following, we provide a generalization of this relation in arbitrary spatial dimensions 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 /4 = π/2ω after the quench, Eq. (36) reduces to W r, p; Using Eq. (9), the spatial density at this time is Making the change of variable u = −p/(mω), the above integral in d dimensions becomes ρ r; and according to Eq. (10), the right-hand side of Eq. (53) essentially represents the initial momentum density, that is 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 /4 after the ballistic expansion in the harmonic trap with frequency ω.

IV. CONCLUDING REMARKS
In this paper we have studied the non-equilibrium dynamics in phase space generated by a sudden change of the Hamiltonian in a quantum system, through the analysis of the Wigner function. For the case of two attractive particles we 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δ (x), a > 0, b ∈ R was considered and we have calculated the two-particle Wigner function. At large values of momentum p, we have found that the associated momentum distribution scales as 1/p 2 , while in the absence of δ interaction, this density scales, in this case of pure zero-range delta interaction, as the well known law 1/p 4 .
We have generalized to arbitrary dimensions d a derivation, by using alternative method, of a relationship shown recently in one dimension between the Wigner function at time t and its initial value following a sudden change of the harmonic trap of noninteracting particles. We have exploited our generalized relation, through the use of inverse Wigner transformation to obtain in d dimensions the scaling law satisfied by the one-body density matrix in real space. Using the generalized relation in Wigner phase space for the considered quench, we have shown that the initial momentum density of a system of particles (interacting or noninteracting) is exactly mapped for a given value of the momentum p =mωr to the spatial density at position r at time a quarter of time period, t = T /4, after the ballistic expansion in the harmonic trap with frequency ω. It should be noted that our method can be easily adapted to deal with physical situations corresponding to quasi-1d or quasi-2d configurations.
An interesting extension of this dynamically situation would be to study the problem of a Lieb-Liniger gas at finite repulsion strength. Work in this direction is in progress.