Quenches from gaussian bosons to Tonks-Girardeau gas: stationary states and effects of a confining potential
QQuenches from bosonic gaussian initial states to the Tonks-Girardeau limit:stationary states and effects of a confining potential
Alvise Bastianello, Mario Collura, and Spyros Sotiriadis
3, 4 SISSA & INFN, via Bonomea 265, 34136 Trieste, Italy The Rudolf Peierls Centre for Theoretical Physics,Oxford University, Oxford, OX1 3NP, United Kingdom. Institut de Math´ematique de Marseille, (I2M) Aix Marseille Universit´e, CNRS,Centrale Marseille, UMR 7373, 39 rue F. Joliot Curie, 13453, Marseille, France University of Roma Tre, Department of Mathematics and Physics, L.go S. L. Murialdo 1, 00146 Roma, Italy
We consider the general problem of quenching an interacting Bose gas from the noninteractingregime to the strongly repulsive limit described by the Tonks-Girardeau gas, with the initial statebeing a gaussian ensemble in terms of the bosons. A generic multi-point correlation function inthe steady state can be fully described in terms of a Fredholm-like determinant suitable both fora numerical and an analytical study in certain limiting cases. Finally, we extend the study tothe presence of a smooth confining potential showing that, in the thermodynamic limit, the timeevolution of the two-point function can be mapped to a classical problem.
I. INTRODUCTION
The recent advances in experimental techniques con-cerning the realm of cold atomic gases made it ur-gent to reach a better understanding of the behavior ofclosed quantum systems when driven out of equilibrium.A paradigmatic out-of-equilibrium protocol consists inthe so called quantum quench , i.e. the system starts inthe ground state, an excited state or a thermal state ofits Hamiltonian and then some parameter of the Hamil-tonian is abruptly changed, resulting in non-trivial dy-namics. Despite the fact that the whole system is closedand so the time evolution is non dissipative, finite sub-systems are expected to relax towards a steady state andimmense efforts have been devoted to investigating thishypothesis .Particularly appealing is the world of one dimensionalsystems, due to the existence of integrable models , i.e.those that possess infinitely many local integrals of mo-tion in addition to the Hamiltonian which allow an ex-act solution of the dynamics . The presence of infinitelymany conserved quantities makes the steady state retainadditional information of the initial state compared tothe non-integrable case so that the Gibbs Ensemble needsto be modified to a Generalized Gibbs Ensemble (GGE)density matrix ρ GGE = 1 Z e − (cid:80) j λ j I j , (1)where I j are the relevant conserved charges and λ j theassociated Lagrange multipliers. The correctness of theGGE prediction has been verified in many different cases and recent works focused in understanding the local-ity features of the charges to be included in the GGEconstruction .Despite several achievements in the study of integrablemodels out-of-equilibrium, apart from free models ana-lytic results have been derived only for restricted classesof observables and initial states and at the price of tremendous efforts . Recent investigations probed inho-mogeneous quench protocols, obtaining analytical resultsin the framework of transport problems .In between the complexity of quenching a truly in-teracting integrable model and the simplicity of a freemodel lies a special interaction quench in the well knownLieb Liniger model . Its Hamiltonian describes a one-dimensional gas of non-relativistic bosons with mass m interacting through a density-density contact potential ofstrength cH LL = (cid:90) dx (cid:20) ∂ x ψ † ( x ) ∂ x ψ ( x )2 m + c ψ † ( x ) ψ † ( x ) ψ ( x ) ψ ( x ) (cid:21) . (2)Above, the scalar fields ψ † ( x ) satisfy the canonicalbosonic commutation relations [ ψ ( x ) , ψ † ( y )] = δ ( x − y ).For arbitrary values of c the model is a truly interact-ing integrable model, however in two special limits thedynamics is that of free particles: at c = 0 the modelreduces to free bosons, while in the c → ∞ limit knownas Tonks-Girardeau gas (TG) the bosons behave as free but impenetrable particles that are equivalent to freefermions by means of a Jordan Wigner transformation.Compared with many other quenches among free theo-ries, the modes at c = 0 and c = ∞ are not connectedthrough a simple linear relation.This quench protocol has already been successfully ad-dressed in the literature, choosing as initial state the BECstate and analytical expressions for the steady state cor-relators were obtained . Subsequently, these results weregeneralized to include a hard wall trap . The presenceof simple analytical results for an experimentally relevantquench protocol is certainly appealing, however the abovementioned works concern a rather peculiar initial state.Among the flaws of the BEC as pre-quench state are: (i) the total absence of any space dependence in the cor-relation functions (except when a confining potential ispresent); (ii) the strong dependence on the trap shape which makes the proper thermodynamic limit difficult todefine. In this sense, the BEC state is rather patholog- a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a y ical and therefore, it is worth to investigate a differentclass of states which hold a more realistic behavior. Onthe other hand, the initial state is still required to besimple enough to allow an analytical treatment. A goodcompromise that meets both requirements is given bythe gaussian (in terms of the bosons) ensembles, whichare the initial states we consider in this paper. In par-ticular, we present a detailed analysis of the multipointcorrelation functions in the steady state attained at latetimes, with emphasis on the large distance decaying prop-erties. We initially consider an homogeneous system andsubsequently a thermodynamically large trap, providingclosed expressions for the correlators suitable both fora numerical and analytical study. The correlation func-tions are notoriously difficult to be obtained in quenchesin truly interacting integrable models: the simplificationof quenching towards the Tonks-Girardeau limit makesthese quantities more easily accessible. While our resultsapply to arbitrary gaussian (in terms of the bosons) ini-tial states, we consider the prototypical example of the free thermal state , i.e. an ensemble with Bose-Einsteindistribution: the density-density correlation function inthe steady state is proven to decay exponentially, the de-cay length being a non trivial function of the temperatureand chemical potential. When the zero temperature limitis taken, the decay behavior changes from exponentialto algebraic: actually, in the zero temperature limit the exact correlation function is computable. This result isstriking different from what is obtained quenching from aBEC , where an exponential decay is observed, but thisis attributed to the fact that the finite temperature GibbsEnsemble is not continuously connected to the BEC stateas long as one dimension is concerned. Our discussion isorganized as follows. In Section II we present the generaltechnology needed to handle this quench protocol: underthe condition of initial clustering of fermionic correla-tors, dephasing among the free fermionic modes of thepost quench dynamics leads to a gaussian steady state in terms of the fermions (i.e. the charges appearing inthe GGE are quadratic in terms of the fermionic opera-tors). Because of the free dynamics connecting the initialstate with the steady state, the open technical problemsreduce to (i) computing the initial fermionic correlatorson a bosonic gaussian ensemble and (ii) computing backthe bosonic correlators on the fermionic gaussian steadystate. We derive closed expressions in terms of Fredholmdeterminants to fulfill this task. These results are em-ployed in Section III to study both numerically and ana-lytically quenches in homogeneous systems, i.e. the par-ticles are confined in a thermodynamically large intervalwith periodic boundary conditions. While being an im-portant benchmark, homogeneous systems with periodicboundary conditions are a crude approximation of real-istic experimental systems, in which an inhomogeneousconfining potential is always present: Section IV extendsour study to the presence of a trap. We show how in thethermodynamic limit (which is necessary for the systemto exhibit relaxation) a semiclassical time evolution natu- rally emerges and allows us to consider traps of arbitraryshapes. Our final conclusions are gathered in Section V,while Appendix A-B contain some technical details insupport to the main text. II. QUENCH SET UP
In the infinite repulsive regime, the Lieb Liniger gas (2)can be described as a free model of hard-core bosons Ψ( x ) H LL = (cid:90) dx ∂ x Ψ † ( x ) ∂ x Ψ( x )2 m , c = + ∞ . (3)The infinite repulsive term prevents the possibility ofhaving two particles in the same position and this canbe achieved through a definition of the hard-core bosonicfields Ψ( x ), obtained projecting the standard bosonicfields ψ ( x ) on the subspace with at most one particlein position x Ψ( x ) = P x ψ ( x ) P x . (4)Above, P x = | x (cid:105)(cid:104) x | + | x (cid:105)(cid:104) x | is the projector on thesubspace with no more than a particle in position x . Theconstraint is evident from the commutation rules satisfiedby the hard-core bosonic fields, i.e. { Ψ( x ) , Ψ( x ) } = 0 , { Ψ( x ) , Ψ † ( x ) } = 1 (5)at the same position and the standard commutation rules[Ψ( x ) , Ψ( y )] = [Ψ( x ) , Ψ † ( y )] = 0 whenever x (cid:54) = y .Actually, a correct definition of the hard-core bosonsrequires a proper lattice regularization (see AppendixA for further details): in order to consider directlythe continuum limit it is convenient to introduce thefermionic field Φ( x ) through a Jordan-Wigner transfor-mation Φ( x ) = exp (cid:0) iπ (cid:82) x dz Ψ † ( z )Ψ( z ) (cid:1) Ψ( x ). While thehard-core bosons need a lattice regularization, on thefermionic fields the continuum limit can be taken andstandard commutation rules { Φ( x ) , Φ † ( y ) } = δ ( x − y ) arerecovered. The strong-interacting Hamiltonian becomingthus quadratic in terms of those fermionic fields.Both models, namely the non-interacting and thestrong-interacting, are therefore related via a non-local transformation which nevertheless makes a well-established connection. Whenever we need to pass fromone theory to the other, e.g. via a quench protocol, theproblem reduces to compute the multipoint correlators ofthe post-quench fields in terms of the pre-quench ones.This dictionary is provided by the underling regularizedhard-core bosons, whose normal ordered correlation func-tions can be computed as if the hard core bosons werecanonical bosonic operators: the details are reported inAppendix A, thereafter we limit ourselves to quote thenecessary results. For example, the fermionic multipointcorrelation function can be expressed in terms of bosonicfields as (cid:104) Φ † ( y ) . . . Φ † ( y l )Φ( x l ) . . . Φ( x ) (cid:105) = (6) (cid:104) : ψ † ( y ) . . . ψ † ( y l ) e − (cid:82) D dz n ( z ) ψ ( x l ) . . . ψ ( x ) : (cid:105) , where, without loss of generality, we assumed y < · · · As first application of the derived expressions we con-sider the quench protocol in the homogeneous case to-wards the Tonks gas. For the seek of concreteness, we of-ten refer to the free thermal ensemble as the pre-quenchstate C ( x, y ) = (cid:90) dk π e ik ( x − y ) e β (cid:16) k m − µ (cid:17) − , (16)where the chemical potential µ is fixed requiring C ( x, x ) = n , with n the particle density. The steady stateinformation is completely encoded in the initial two pointfermionic correlator, that is conserved through the timeevolution. While the numerical evaluation of Eq. (14) canbe easily implemented, its analytical solution can be ob-tained only in very special cases, since the presence ofthe finite domain D prevents the possibility of diagonal-izing the operators in the Fourier space. However, Eq.(14) is still feasible of proper analytical approximationsin several regimes. a. Short distances.— In the short distance regimeor in the small density limit, a satisfactory approximationof Eq. (14) can be obtained from the series expansion interms of C D . For example, assuming low density and thatthe two point bosonic correlators decays fast enough toensure the convergence of the integrals in Eq. (10-11) forinfinite domain D , we can readily write (cid:104) Φ † ( y )Φ( x ) (cid:105) = C ( y, x ) e − | x − y | ( n + O ( n ) ) + O ( n ) . (17) b. Small temperatures.— When the two point cor-relator is approximately constant C ( x, y ) (cid:39) C ( x, x ) = n over the whole domain D . This is an important case,since this condition is met in the zero temperature limit,at fixed density, of the free thermal ensemble: in thelow temperature limit C ( x, y ) (cid:39) ne − α | x − y | with α = m ( nβ ) − . Thus in the small temperature limit α ap-proaches zero and C is approximatively constant on anincreasing range of distances. In this limit the followingtwo point fermionic correlator is immediately derived (cid:104) Φ † ( y )Φ( x ) (cid:105) = n (1 + 2 n | x − y | ) . (18)This correlator greatly differs from the findings in Ref.14 for the BEC case, where the decaying has been foundto be purely exponential rather than algebraic. The tworesults are not in contradiction, since the condensationof one dimensional systems is realized only exactly atzero temperature, thus the zero temperature limit of thethermal ensemble does not reproduce the BEC state. c. Large distances.— When the distance of the twofermionic operators in the two point correlator far ex-ceeds the typical decay length of the bosonic correlator C ,we can extract the asymptotics of the two point fermioniccorrelator. We start considering the string contributiondet D (1 + 2 C D ) − = e − (cid:82) yx dz [log(1+2 C D )]( z,z ) , (19)where we already specialized the domain D = ( x, y ).In the limit of a very large interval, most of the con-tribution to the integral comes from coordinates far fromthe edges of the domain. Thus, in first approximation,we can compute [log(1 + 2 C D )]( z, z ) as if the domain D were the whole real axis and diagonalize the operator inthe Fourier space. This approximation gives the extensivepart of the integraldet D (1 + 2 C D ) − = e −| x − y | (cid:82) dk π log [ C ( k ) ] + O ( | x − y | ) , (20)where ˜ C ( k ) is the Fourier transform of the two pointbosonic correlator. The same approximation can be jus-tified on the whole two point fermionic correlator, beingreliable only in the asymptotic regime (cid:104) Φ † ( y )Φ( x ) (cid:105) ∝ e −| x − y | (cid:82) dk π log [ C ( k ) ] (21) × (cid:90) dk π e ik ( x − y ) ˜ C ( k )1 + 2 ˜ C ( k ) . In the case of the thermal ensemble, we have (cid:104) Φ † ( y )Φ( x ) (cid:105) ∝ e − ξ | x − y | (cid:18) e − λ | x − y | λ + c.c. (cid:19) , (22)where the parameters ξ and λ are ξ = (cid:90) dk π log coth (cid:20) β (cid:18) k m − µ (cid:19)(cid:21) , (23) λ = √ m (cid:118)(cid:117)(cid:117)(cid:116) | µ | + (cid:115) µ + π β + i (cid:118)(cid:117)(cid:117)(cid:116) −| µ | + (cid:115) µ + π β . The analysis of the asymptotes displays an oscillating be-havior due to the imaginary part of λ , however in view x -14 𝛽 = 0.5 𝛽 = 1 𝛽 = 2 𝛽 = 8 𝛽 = 64 𝛽 = ∞ z 𝛷 † ( z ) 𝛷 ( ) | ⟨ ⟩ | m = 1, n = 1 𝛽 = 0.5 𝛽 = 1 𝛽 = 2 𝛽 = 8 𝛽 = 64 𝛽 = ∞ k n ( k ) m = 1, n = 1 ~ x -6 𝛽 = 0.5 𝛽 = 1 𝛽 = 2 𝛽 = 8 𝛽 = 64 𝛽 = ∞ z 𝜓 † ( z ) 𝜓 ( ) m = 1, n = 1 𝛽 = 0.5 𝛽 = 1 𝛽 = 2 𝛽 = 8 𝛽 = ∞ FIG. 2: (Left) Stationary fermionic two-point correlation function after an interaction quench ( c = 0 → c = ∞ ) of a Bose gasinitially prepared in a thermal state with inverse temperature β . Symbols are the numerical evaluation of Eq. (14) with initialbosonic correlator given by (16). In the limit of zero temperature, the correlation function approaches the analytic result (18) (fullred line). However, for any finite temperature, the large-distance exponential decay is captured by the asymptotic expansion (22)(dashed lines). (Center) Stationary momentum density distribution of the fermions obtained by Fourier transform (cid:104) Φ † ( z )Φ(0) (cid:105) .(Right) Stationary bosonic two-point correlation function obtained by numerical evaluation of Eq. (15) with initial fermioniccorrelator given by (cid:104) Φ † ( z )Φ(0) (cid:105) . Notice the large distance exponential decay. Nevertheless, such decay is not exponential for alldistances. Indeed, by comparing the correlator with the exact exponential decay which matches the short distance expansion,namely n e − n | x | (dot-dashed black line), one reveals an appreciable discrepancy (see the inset for a zoom on the short distances). of our crude approximations only the overall exponentialdecay (cid:104) Φ † ( y )Φ( x ) (cid:105) ∝ e − ( ξ + (cid:60) ( λ )) | x − y | is expected to be re-liable. In Fig. 2 we plot the numerical computation of thestationary two point fermionic correlator after a quenchfrom a thermal bosonic ensemble (16) and compare itwith our analytical predictions. As already stated before,we only report the pure exponential decay. Notwithstand-ing, we checked the predicted period 2 π/ (cid:61) ( λ ) againstthe numerical results and, quite remarkably, we founda very good agreement. However, the asymptotic expan-sion does not give a correct prediction for the oscillations’phase.We also analyze the momentum density distribution˜ n ( k ) = (cid:82) dz e − ikz (cid:104) Φ † ( z )Φ(0) (cid:105) , which turns out to be verydifferent from the fermionic momentum distribution atthermal equilibrium. In general, for any finite tempera-ture, the momentum density distribution is an analyticfunction for all momenta. Only in the zero temperaturelimit, the algebraic decay of the fermionic two-point func-tion affects the short-distance behavior of the momen-tum distribution which reads ˜ n ( k ) (cid:39) − π | k | / n (for β → ∞ ). On the contrary, the large momentum behav-ior ˜ n ( k ) (cid:39) n /k it turns out to be independent of thetemperature, being fixed by the non analyticity of thefermionic correlator (cid:104) Φ † ( z )Φ(0) (cid:105) (cid:39) n − n | z | in z = 0.Let us point out that such non-analytic behavior in thefermionic two point function is somehow universal sinceit only depends on the properties of the initial bosoniccorrelator in the vicinity of z = 0, namely C ( x, x ) = n and ∂ x C ( x, y ) | y = x = − ∂ y C ( x, y ) | y = x = 0.Finally, using the steady-state of the two-pointfermionic correlation function, we can numerically eval-uate the stationary bosonic correlators (cid:104) ψ † ( z ) ψ (0) (cid:105) via Eq. (15). The Bose gas in the steady state is fully char-acterized by a two-point function which exhibits a nearlyexact exponential decay for all temperatures, even for β = ∞ , even tough the stationary fermionic two-pointfunction decays algebraically. Actually, the exponentialdecay of the bosonic correlator in the steady state can beargued by mean of the symmetric analysis performed toextract the asymptotic decay of the two-point fermioniccorrelator. Interestingly, for all temperatures, we can eas-ily extract the short distance behaviour of the bosonictwo-point function which, in particular, coincides withthe fermionic one, i.e. (cid:104) ψ † ( z ) ψ (0) (cid:105) (cid:39) (cid:104) Φ † (0)Φ(0) (cid:105) + ∂ z (cid:104) Φ † ( z )Φ(0) (cid:105)| z = n − n | z | .Lastly, we remark that even though the asymptoticbehavior of the two point fermionic correlator does notsuffice to determine the decay of an arbitrary bosoniccorrelator because of the presence of the string, this isinstead possible when the string is absent. This is thecase, for example, in the density-density correlator, beingthe density equivalently expressed in terms of either thefermions or the bosons. (cid:104) n ( x, t ) n ( y, t ) (cid:105) steady state = n − |(cid:104) Φ † ( x )Φ( y ) (cid:105)| (24)Thus, the density-density correlator (connected part) de-cays twice as fast as the two point fermionic correlator(22) and it can be computed exactly in the zero temper-ature limit (18). 𝛕 = 0 𝛕 = 𝝅 /4 𝛕 = 𝝅 /2 𝛕 = 3 𝝅 /4 𝛕 = 𝝅𝛕 = 0 𝛕 = 2 𝝅 𝛕 = 4 𝝅 𝛕 = 16 𝝅 𝛕 = 32 𝝅 v ( x ) = x v ( x ) = x FIG. 3: Color density plot of the Wigner distribution F τ ( p, q ) vs ( q, p ) ∈ [ − , × [ − , τ and tworepresentative traps (here m = 1 and n = 1). The system is initially prepared at zero temperature with F ( p, q ) given byEq. (33). In the harmonic case each phase-space trajectory is characterized by the same period thus preserving the shape ofthe initial Wigner distribution whose support simply rotates in time. In the anharmonic case, the typical dynamical period T ( p ) (cid:39) (cid:82) q ∗ dq/ (cid:112) p − mV ( q ) (with q ∗ being the de real positive root of V ( q ∗ ) = p / m ) depends on the initial energy, andthe δ -support of the Wigner distribution folds over time. For large times, it tends to cover the entire phase plane. IV. QUENCHING IN A THERMODYNAMICTRAP The expression (14) for the correlators does not requiretranslational invariance, but only gaussianity of the ini-tial state. Thus, we can generalize the previous calcula-tions to the more realistic situation in which a confiningpotential is present and replace the Hamiltonian (2) with H LL → H LL + (cid:90) dx V ( x ) ψ † ( x ) ψ ( x ) , (25)where V ( x ) is a trapping potential. Differently from thehomogeneous case, the two point fermionic correlator ex-hibits now a non-trivial time evolution.In presence of the trap, the discrete spectrum preventsthe relaxation to a steady state, that can be recoveredonly in a proper thermodynamic limit: for this reason, weintroduce a thermodynamic length-scale L . Parametriz-ing the confining potential as V ( x ) = v ( xL − ), where v is assumed to be a smooth function kept fixed in thethermodynamic limit, the eigenvalues approach a con-tinuum distribution when L → ∞ . While consideringthe limit of large trap, we must also rescale the num-ber of particles consistently and keep the total density n = L − (cid:82) dz C ( z, z ) constant.The thermodynamic limit hugely simplifies the sub-sequent calculations and permits to consider a generictrap. First of all, in the thermodynamic limit the ther-mal initial state is described by the Local DensityApproximation (LDA). In fact, the two point correlatordecays on a typical scale (finite for any non zero tempera-ture) much smaller than the thermodynamic length char-acterizing the variation of V . It is convenient to change variables in the two point bosonic correlator C ( y, x ) → C ( s, q ), defining s = y − x and q = ( x + y ) / L . In theLDA we have C ( s, q ) = (cid:90) dk π e iks e β (cid:16) k m + v ( q ) − µ (cid:17) − . (26)As is clear from the general expression (14), the LDAon the bosonic correlators implies LDA for the fermioniccorrelators as well, i.e. F ( x, y ) → F ( s, q ).Following the time evolution after the quench, we canclearly distinguish two time scales.1) On a non-thermodynamic time scale, the systembehaves as if it was locally homogeneous, therefore it ex-hibits gaussification of the correlators on the same timescale as the homogeneous system. In this first phase, thetwo point correlator does not evolve: even though the en-semble is gaussian, the steady state has not been reachedyet.2) Relaxation can be realized only on a time scale suf-ficient for the system to explore the whole trap, thustimes t ∝ L must be considered. The two point correla-tor evolves under the Schr¨odinger equation, that in termsof the rescaled variables reads i∂ t F t = (cid:26) ∂ s ∂ q Lm − (cid:104) v (cid:16) q + s L (cid:17) − v (cid:16) q − s L (cid:17)(cid:105)(cid:27) F t . (27)For any finite time t , in the limit L → ∞ the twopoint correlator does not evolve. However, rewriting thedifferential equation in terms of the rescaled time τ = t/L and Taylor expanding the potential v , we get a non trivialtime evolution i∂ τ F τ = (cid:20) ∂ s ∂ q m − sv (cid:48) ( q ) (cid:21) F τ (28)where we assumed τ (cid:28) L (i.e. t (cid:28) L ). The above equa-tion is associated with a classical evolution for the twopoint correlator. This is unveiled in terms of the wellknown Wigner distribution , that amounts to a partialFourier transform of the correlator˜ F τ ( p, q ) = (cid:90) ds e − ips F τ ( s, q ) . (29)In this new notation, Eq. (28) can be written in the formof a classical Liouville equation in terms of q and its con-jugated momentum p∂ τ ˜ F τ + { H cl , ˜ F τ } p,q = 0 , (30)where { , } p,q are the classical Poisson brackets and H cl isthe classical Hamiltonian H cl = p / m + v ( q ). In fact, thethermodynamic limit we are considering is equivalent toa semiclassical limit on the time evolution of the Wignerdistribution .If the frequency of the solution of the classical equa-tion of motion does not have a trivial dependence on theenergy, the correlator in the coordinate space reaches asteady state thanks to the dephasing between the differ-ent energy shells. The steady state is of course describedin terms of the energy density of the initial Wigner dis-tributionlim τ →∞ F τ ( s, q ) = (cid:90) dp π e ips (cid:15) (cid:18) p m + v ( q ) (cid:19) , (31)where the energy density is defined as (cid:15) ( E ) = (cid:82) dpdq δ (cid:16) p m + v ( q ) − E (cid:17) ˜ F ( p, q ) (cid:82) dpdq δ (cid:16) p m + v ( q ) − E (cid:17) . (32)Based on the integrability of the free fermionic model,the steady state is expected to be described by a GGEconstructed out of the modes of the trap: in AppendixB we show that the quantum mechanical prediction co-incides with the above semiclassical result, as it shouldbe. As already mentioned, relaxation is caused by the de-phasing between the different classical trajectories. Whilethis is the common scenario for a generic trap, this is nolonger true in the parabolic case: being the period of theoscillations independent from the initial conditions, de-phasing is absent and the two point correlator exhibitsrevivals. Amusingly, notice that for a parabolic potentialEq. (27) coincides with Eq. (28) without taking the ther-modynamic limit, thus the Wigner distribution evolvesclassically even for finite harmonic traps.As an example, we explicitly analyze the zero temper-ature limit. In this limit the LDA of the initial fermionictwo-point correlation function reduces to F ( s, q ) = δ ( q ) n (1 + 2 n | s | ) , (33)and the Wigner distribution evolves as˜ F τ ( p, q ) = (cid:90) dp δ ( q − q τ ( p )) δ ( p − p τ ( p )) ˜ f ( p ) , (34) where ˜ f ( p ) = (cid:82) ds e − ip s n/ (1 + 2 n | s | ) , and( q τ ( p ) , p τ ( p )) is the solution of the classical equationof motion with initial condition (0 , p ). Therefore,exploiting the normalization of the Dirac delta function,we can integrate over the final momentum p , obtainingfor the correlation function F τ ( s, q ) = (cid:90) dp π δ ( q − q τ ( p ))e isp τ ( p ) ˜ f ( p ) (35)= (cid:88) j π e isp τ ( p ) ˜ f ( p ) | ∂ p q τ ( p ) | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p = p j ( q,τ ) , where the sum runs over all the real roots p j ( q, τ ) ofthe equation q τ ( p j ) = q . In Fig. 3 we plot the Wignerdistribution at different times τ , and for two represen-tative confining potentials, namely v ( x ) = x / v ( x ) = x / m = 1and n = 1). In the harmonic case, the Wigner distribu-tion exhibits a trivial periodic dynamics and the particledensity never reaches a stationary profile, continuouslyoscillating between δ ( q ) and ˜ f ( q ) / π .On the contrary, in the anharmonic case, the number ofroots of the equation q τ ( p j ) = q increases with time andthe root density approaches a continuous distribution;the system therefore relaxes to the stationary state whosecorrelation function is given by Eq. (31).In Fig 4 we compare the steady-state predictions withthe correlation function evaluated numerically using Eq.(35), in the presence of the anharmonic confining po-tential v ( x ) = x / 4. In particular, we plot the particledensity versus q as well as the two-point correlator as afunction of s = y − x for fixed q = 0 and q = 2. As ex-pected, for large time the numerical data approach thestationary profiles. V. CONCLUSION In this work we studied the equilibration mechanismof 1D Bose gases in the Tonks-Girardeau limit startingwith initial states that are gaussian in the bosonic fields,in the homogeneous case as well as the general inhomo-geneous case of an arbitrary confining trap. Despite thenoninteracting (in terms of the fermions) dynamics, thisproblem is nontrivial due to the nonlocal nature of theJordan-Wigner transformation which links the bosonicand fermionic fields. By solving the problem of extractingfermionic correlations from the bosonic ones in a Gaus-sian state and vice versa, we derived exact closed formu-las for the steady state correlations and analysed themboth numerically and analytically in several asymptoticlimits.When a confining potential is considered, the relax-ation of local observables can be attained only in a properthermodynamic limit, where the local density approxi-mation is applicable. In this limit, the evolution is sepa-rated on two time scales: initially the system behaves as -4 -2 0 2 400.10.20.30.40.5 𝛕 = 2 𝝅𝛕 = 4 𝝅𝛕 = 16 𝝅𝛕 = 32 𝝅 v ( x ) = 14 x q F 𝛕 ( , q ) m = 1, n = 1 -4 -2 0 2 400.10.20.30.4 𝛕 = 𝝅 /2 𝛕 = 2 𝝅𝛕 = 4 𝝅𝛕 = 16 𝝅𝛕 = 32 𝝅 v ( x ) = 14 x s | F 𝛕 ( s , ) | m = 1, n = 1 -2 -1 0 1 200.050.10.15 𝛕 = 𝝅 /2 𝛕 = 2 𝝅𝛕 = 4 𝝅𝛕 = 16 𝝅𝛕 = 32 𝝅 v ( x ) = 14 x s | F 𝛕 ( s , ) | m = 1, n = 1 FIG. 4: Zero temperature two-point correlation function for different times after an interaction quench in the anharmonic trap.In the large time limit the correlator approaches the steady-state prediction of Eq. (31) (black dashed lines). if it were homogeneous with a subsequent gaussianifica-tion. Thereafter, the particles explore the whole trap ona thermodynamic time scale with an emergent semiclas-sical evolution for the two point fermionic correlator, interms of which all the expectation values of local observ-ables are determined. VI. ACKNOWLEDGMENTS This work was supported by the European Union’sHorizon 2020 research and innovation program under theMarie Sklodowska-Curie Grant Agreement No. 701221(M.C.). SS acknowledges support from the A*MIDEXproject Hypathie (no. ANR-11-IDEX- 0001-02) and fromthe LIA - LYSM (Laboratoire Ypatia des SciencesMath´ematiques). Appendix A: Correlators of fermions on the initialstate The key point in the analysis of our quench protocol isthe mapping between bosonic and fermionic correlators,achieved by mean of the Jordan Wigner transformation(6). This mapping passes through the hard core bosons,whose rigorous treatment requires a lattice regulariza-tion. This appendix is devoted to this issue, in particularwe show that even though fermionic fields are definedin terms of the hard core bosons, their correlation func-tions can be computed as if the hard core bosons werestandard bosonic fields (after the expressions have beennormal ordered). In the case when the initial state is cho-sen to be the Bose-Einstein condensate the explicit checkhas already been performed in Ref. 14, however the sameconclusions hold even for more general states. Througha simple generalization of the following, it is possible toshow the inverse relation, i.e. that bosonic correlators areexpressed in terms of fermionic correlators exchangingthe role of bosons and fermions in Eq. (6). Here we test expectation values on pure states witha fixed number of particles, but the generalization overmore general density matrices (as the thermal ensemble)is trivial. Consider then a generic state | S (cid:105) with n parti-cles | S (cid:105) = (cid:90) d n x W n ( x , ..., x n ) ψ † ( x ) ...ψ † ( x n ) | (cid:105) , (A1)where W n is the symmetric (normalizable) wavefunction,that we will assume to be smooth. We introduce a latticespace ∆ and discretize the space x → ∆ j , the continuousbosonic fields are replaced with the lattice counterpart ψ (∆ j ) → ∆ − / ψ j | S ∆ (cid:105) = ∆ n/ (cid:88) j ,...,j n W n (∆ j , ..., ∆ j n ) ψ † j ...ψ † j n | (cid:105) , (A2)where the lattice regularized fields follow standard com-mutation rules [ ψ j , ψ † j (cid:48) ] = δ j,j (cid:48) . To be precise, the latticewavefunction differs from the continuous one for terms O (∆), but we drop them since they are inessential inthe limit in which we are interested in. The hard corebosons are obtained from the standard bosons by meanof a projector that excludes the presence of more thanone particle at each lattice siteΨ j = P j ψ j P j , (A3)where P j = | j (cid:105)(cid:104) j | + | j (cid:105)(cid:104) j | . By mean of a straightfor-ward application of this definition, the hardcore bosonicfields are immediately shown to satisfy the following com-mutation rules { Ψ j , Ψ † j } = 0 , Ψ j = 0 (A4)[Ψ j , Ψ j (cid:48) ] = [Ψ j , Ψ † j (cid:48) ] = 0 , j (cid:54) = j (cid:48) (A5)With the hard core bosons we define the fermions onthe lattice through the Jordan Wigner transformationΦ j = exp iπ (cid:88) j (cid:48) In Section IV we considered the quench protocol ina thermodynamic trap and we showed the emergenceof a classical picture for the Wigner distribution, withthe subsequent classical relaxation to a steady state. Onthe other hand, based on the integrability of the model,we expect relaxation to a GGE constructed out of themodes of the trap. This appendix is devoted to show-ing the equivalence of the two descriptions, i.e. that theGGE predictions reduce, in the thermodynamical limit,to the classical steady state (31-32). This consistencycheck can be achieved thanks to the fact that the WKBapproximation becomes exact in the TDL limit we areinvestigating.We start presenting a brief review of the WKB methodsufficient for our purposes, then use it to compute theGGE predictions in the thermodynamic limit. 1. WKB description of the thermodynamic trap Consider the one particle Shroedinger equation in onedimension: − ∂ x ϕ ( x )2 m + V ( x ) ϕ ( x ) = Eϕ ( x ) . (B1)The trapping potential is assumed in the form of SectionIV, i.e. V ( x ) = v ( xL − ). For simplicity, we assume v tobe symmetric and with an unique minimum in x = 0.In the limit of small variation of the potential the wave-function can be approximated through the WKB solu-tion. Let ξ n > E n = v ( ξ n ) ( E n is the energy of the n th energy level), then the wavefunction oscillates inside the0trap ( | x | < Lξ n ) ϕ n ( x ) = N n cos (cid:16) πn + (cid:82) x dy (cid:112) m ( E n − v ( y/L )) (cid:17) [2 m ( E n − v ( x/L ))] / (B2)and is exponentially damped outside of it. For example,on the right hand side x > Lξ n we have ϕ n ( x ) ∝ m ( v ( x/L ) − E n )] / e − (cid:82) xLξn dy √ m ( v ( y/L ) − E n ) . (B3)In Eq. (B2) N n is a normalization constant to ensurethe eigenfunction to be normalized to unity. The energylevels satisfy the WKB quantization condition (cid:90) ξ n − ξ n dq (cid:112) m ( E n − v ( q )) = πL (cid:18) n + 12 (cid:19) . (B4)As we already said, the WKB approximation is reliablewhen the variation of the potential is much smaller thanthe energy difference between E and the potential, morequantitatively this amounts to ask (cid:12)(cid:12)(cid:12)(cid:12) ∂ x V ( E n − V ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:28) , (cid:12)(cid:12)(cid:12)(cid:12) ( ∂ x V ) ( E n − V ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:28) . (B5)Replacing V with the rescaled potential v (and using therescaled variable q = xL − ) we see that the WKB ap-proximation becomes exact in the TDL1 L (cid:12)(cid:12)(cid:12)(cid:12) v (cid:48)(cid:48) ( q )( E n − v ( q )) (cid:12)(cid:12)(cid:12)(cid:12) (cid:28) , L (cid:12)(cid:12)(cid:12)(cid:12) ( v (cid:48) ( q )) ( E n − v ( q )) (cid:12)(cid:12)(cid:12)(cid:12) (cid:28) , (B6)with the exception of the turning points, i.e. the zeros ofthe denominators of the above expressions. However, call-ing δq the distance from the turning point, the region ofvalidity of the WKB approximation is readily estimatedas 1 | v (cid:48) ( ξ n ) | L (cid:28) | δq | . (B7)Thus, in the rescaled variable, the region in which theWKB approximation fails is shrunk to zero as δq ∝ L − / . Some care must be used stating that the WKB be-comes exact in the thermodynamic limit, since in termsof the original variable x = Lq the region where WKBfails does not vanish, but rather diverges as L / . Indeed,even though the single WKB eigenfunction is not reliablein a neighborhood of the turning point of increasing size,the non extensive scaling of the latter permits, ipso facto,a direct use of the WKB eigenfunctions in several expres-sions.Before passing to analyze the GGE, we need to prop-erly normalize to unity the wavefunction: this task ishugely simplified in the thermodynamic limit, since only the WKB solution inside the trap is needed. To un-derstand this point consider the solution outside of thetrap, let’s choose the right side for simplicity. Using therescaled coordinate q = L − x the wavefunction is expo-nentially damped as e − L (cid:82) qξn dq (cid:48) √ m ( v ( q (cid:48) ) − E ) . (B8)In the limit L → ∞ the above exponential vanishes,meaning that the region where the wavefunction is notzero outside of the trap does not have a thermodynamicextent. Using this fact (and the exactness of the WKBapproximation in the rescaled coordinate), the normal-ization constant N n in (B2) is readily written as N n = (cid:34) L (cid:90) ξ n − ξ n dq m ( E n − v ( q ))] / (cid:35) − = 2 √ m (cid:112) Lρ ( E n ) , (B9)where we neglected subextensive contributions to N − n , ρ ( E ) is simply the volume of the energy shell in the clas-sical phase space ρ ( E ) = (cid:90) dτ dp δ (cid:18) p m + v ( τ ) − E (cid:19) . (B10) 2. The excitation density in the thermodynamiclimit Armed with the WKB approximation, we can startconstructing the thermodynamic limit of the GGE asso-ciated with the trap. As first step, we extract the excita-tion energy from the two point fermionic correlator in theassumption the latter satisfied the LDA approximation,i.e. (cid:104) Φ † ( x )Φ( y ) (cid:105) = F ( x, y ) → F ( s, q ) where s = x − y and q = ( x + y ) / L . F does not contain any further depen-dence on the thermodynamic length L and decays fastenough for large s to justify the forthcoming approxima-tions.Denoting as η † n , η n the creation/annihilation operatorsof the modes of the trap, the number of particles in theenergy level (cid:104) η † n η n (cid:105) can be readily extracted from the twopoint fermionic correlator and the eigenfunctions (cid:104) η † n η n (cid:105) = (cid:90) dxdy ϕ ∗ n ( x ) ϕ n ( y ) F ( x, y ) . (B11)At the price of neglecting corrections vanishing in thethermodynamic limit, we can replace in the above theWKB eigenfunctions and the LDA approximation for thecorrelator. Using (B2) with the truncation of the wave-function outside of the trap, through simple trigonomet-ric identities we get1 (cid:104) η † n η n (cid:105) = 2 mLρ ( E n ) (cid:90) Lξ n − Lξ n dxdy cos (cid:16)(cid:82) xy dx (cid:48) (cid:112) m ( E n − v ( x (cid:48) /L )) (cid:17) [2 m ( E n − v ( x/L ))] [2 m ( E n − v ( y/L ))] F (cid:18) x − y, x + y L (cid:19) (B12)+ 2 mLρ ( E n ) (cid:90) Lξ n − Lξ n dxdy cos (cid:16) πn + (cid:82) x dx (cid:48) (cid:112) m ( E n − v ( x (cid:48) /L )) + (cid:82) y dy (cid:48) (cid:112) m ( E n − v ( y (cid:48) /L )) (cid:17) [2 m ( E n − v ( x/L ))] [2 m ( E n − v ( y/L ))] F (cid:18) x − y, x + y L (cid:19) . Thanks to the LDA approximation for the fermionic cor-relator and under the assumption that decays in the rel-ative distance x − y , the above can be further simplified.In particular, the second term gives contributions sup-pressed in the thermodynamic limit and the first line canbe approximated to (cid:104) η † n η n (cid:105) = (B13) (cid:90) ξ n − ξ n dq (cid:90) ∞−∞ ds cos (cid:16) s (cid:112) m ( E n − v ( q )) (cid:17) [2 m ( E n − v ( q ))] mF ( s, q ) ρ ( E n ) . The above is clearly a smooth function of the energy (cid:104) η † n η n (cid:105) = (cid:15) ( E n ), where (cid:15) ( E ) = 1 ρ ( E ) (cid:90) dpdq δ (cid:18) p m + v ( q ) − E (cid:19) ˜ F ( p, q )(B14)and ˜ F is defined in Eq. (29). Notice that (cid:15) coincides withthe classical expression (32). 3. GGE predictions in the thermodynamic limit The GGE density matrix of the free fermion gas is de-fined over the modes of the trap ρ GGE = e − (cid:80) n λ n η † n η n / Z ,in particular this ensemble is gaussian and it is com-pletely determined by the two point correlator (cid:104) Φ † ( x )Φ( y ) (cid:105) GGE = (cid:88) n ϕ ∗ n ( x ) ϕ n ( y ) (cid:104) η † n η n (cid:105) GGE , (B15)where the GGE expectation value of the occupancy η † n η n is equal to the expectation value on the initial state, that reduces to (B14) in the TDL. Computing (B15) inthe thermodynamic limit simply amounts to replace theeigenfunctions with the WKB expressions and the modeoccupation with (B14). In the TDL limit the energy lev-els become dense, thus we would like to replace the sumover the discrete levels with an integral over the energy,being the density of levels at fixed energy extracted fromthe quantization rule (B4). However, the energy eigen-functions are not smooth functions of the energy, sincethey are even (odd) if the quantum number is even (odd),as it is clear from (B2). Therefore, we first split the sumover the even and odd eigenfunctions, then take the con-tinuum approximation of each one. Doing such an oper-ation and recollecting the terms in an unique integral wereadily find (cid:104) Φ † ( x )Φ( y ) (cid:105) GGE = (B16) (cid:90) dE π m cos (cid:16)(cid:82) xy dx (cid:48) (cid:112) m ( E − v ( x (cid:48) /L )) (cid:17) [2 m ( E − v ( x/L ))] [2 m ( E − v ( y/L ))] (cid:15) ( E ) . Since (cid:15) ( E ) is a smooth function, the above integral decayswhen x and y are pulled far apart and in the limit of large L we can safely approximate (cid:104) Φ † ( x )Φ( y ) (cid:105) GGE = (cid:90) dp π e ips (cid:15) (cid:18) p m + v ( q ) (cid:19) , (B17)where, as usual, s = x − y and q = ( x + y ) / L . Theabove expression matches with the classical computation(31-32), as it should be. M. Greiner et al , Nature , 51-54 (2002); T. Kinoshita,T. Wenger, and D. S. Weiss, Nature , 900 (2006);S. Hofferberth, I. Lesanovsky et al , Nature , 324-327(2007); L. Hackermuller, U. Schneider et al , Science ,1621 (2010); S. Trotzky, Y.-A. Chen et al , Nature Phys. ,325 (2012); M. Gring, M. Kuhnert et al , Science , 1318(2012): M. Cheneau, P. Barmettler et al , Nature , 484(2012); T. Langen, R. Geiger et al , Nature Physics , 640(2013); F. Meinert, M.J. Mark et al , Phys. Rev. Lett. ,053003 (2013); J.P. Ronzheimer, M. Schreiber et al , Phys.Rev. Lett. , 205301 (2013). P. Calabrese and J. Cardy, Phys. Rev. Lett. , 136801(2006); J. Stat. Mech. (2007) P06008. A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat-tore, Rev. Mod. Phys. , 863 (2011). V.E. Korepin, N.M. Bogoliubov, A,G, Izergin, Quan-tum Inverse Scattering Method and Correlation Functions (University Press, Cambridge, 1993). F.H.L. Essler and M. Fagotti, J. Stat. Mech. (2016)064002, Special issue on Quantum Integrability in Out ofEquilibrium Systems , editors P. Calabrese, F.H.L. Esslerand G. Mussardo, J. Stat. Mech. (2016) 064001;M. A. Cazalilla, Phys. Rev. Lett. , 156403 (2006); M. Cramer,C. M. Dawson, J. Eisert, and T. J. Osborne, Phys. Rev.Lett. , 030602 (2008); M. Cramer, J. Eisert, NewJ. Phys. , 055020 (2010); A. Silva, Phys. Rev. Lett. , 120603 (2008); M. Collura, S. Sotiriadis, and P. Cal-abrese, Phys. Rev. Lett. , 245301 (2013); J. Stat. Mech.(2013) P09025; G. Mussardo, Phys. Rev. Lett. , 100401(2013); G. Goldstein and N. Andrei, Phys. Rev. A ,043625 (2014); M. Collura and D. Karevski, Phys. Rev.B. , 214308 (2014); M. Collura, P. Calabrese, F. H. L.Essler, Phys. Rev. B , 125131 (2015); B. Bertini, L.Piroli, P. Calabrese, J. Stat. Mech. (2016) 063102; A. Bas-tianello, A. De Luca, G. Mussardo, J. Stat. Mech. (2016)123104. M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys.Rev. Lett. , 050405 (2007). J.-S. Caux and F.H.L. Essler, Phys. Rev. Lett. , 257203(2013); B. Pozsgay, J. Stat. Mech. P10045 (2014); P. Cal-abrese, F. H. L. Essler and M. Fagotti, Phys. Rev. Lett. , 227203 (2011); M. Fagotti and F.H.L. Essler, Phys.Rev. B , 245107 (2013); F. H. L. Essler, S. Evangelisti,and M. Fagotti, Phys. Rev. Lett. , 247206 (2012). S. Sotiriadis, P. Calabrese, J. Stat. Mech. (2014) P07024; S.Sotiriadis, G. Martelloni, J. Phys. A , 1020 (2016); S. Sotiriadis,arXiv:1612.00373 B. Pozsgay, M. Mesty´an et al , Phys. Rev. Lett. , 117203(2014); M. Mestyan, B. Pozsgay, G. T´akacs, M.A. Werner,J. Stat. Mech. P04001 (2015); E. Ilievski, J. De Nardis,B. Wouters, J-S Caux, F. H. L. Essler, T. Prosen Phys.Rev. Lett. , 157201 (2015); E. Ilievski, M. Medenjak,and T. Prosen, Phys. Rev. Lett. , 120601 (2015); E.Ilievski, M. Medenjak et al , J. Stat. Mech. 064008 (2016);E. Ilievski, E. Quinn, J. De Nardis, M. Brockmann J. Stat.Mech. P063101 (2016); L. Piroli, E. Vernier, J. Stat. Mech.P053106 (2016); S. Sotiriadis, Phys. Rev. A , 031605(2016); F. H. L. Essler, G. Mussardo, and M. Panfil, Phys.Rev. A , 051602(R) (2015); A. Bastianello, S. Sotiriadis,J. Stat. Mech. 023105 (2017) ; F. H. L. Essler, G. Mussardo,M. Panfil, J. Stat. Mech. 013103 (2017); E. Vernier and A.C. Cubero, J. Stat. Mech. 023101 (2017). B. Pozsgay, J. Stat. Mech. P07003 (2013);M. Fagotti andF. H. L. Essler, J. Stat. Mech. P07012 (2013); W. Liu andN. Andrei, Phys. Rev. Lett. , 257204 (2014); J. DeNardis, B. Wouters, M. Brockmann, and J.-S. Caux, Phys.Rev. A , 033601 (2014); M. Fagotti, M. Collura et al , Phys. Rev. B , 125101 (2014); B. Wouters, J. De Nardis et al , Phys. Rev. Lett. , 117202 (2014); L. Piroli, E.Vernier, P. Calabrese, Phys. Rev. B , 054313 (2016). B. Bertini, M. Collura, J. De Nardis, M. Fagotti, Phys.Rev. Lett. , 207201 (2016); O. A. Castro-Alvaredo,B. Doyon, T. Yoshimura, Phys. Rev. X , 041065 (2016);A. De Luca, M. Collura, J. De Nardis, arXiv:1612.07265(2017). E. H. Lieb and W. Liniger, Phys. Rev. , 1605 (1963),J.-S. Caux and R. M. Konik, Phys. Rev. Lett. , 175301(2012); M. Kormos, A. Shashi, Y.-Z. Chou, J.-S. Caux,and A. Imambekov, Phys. Rev. B , 205131 (2013). J. DeNardis, B. Wouters, M. Brockmann, and J.-S. Caux, Phys.Rev. A , 033601 (2014); J. De Nardis, J.-S. Caux, J.Stat. Mech. P12012 (2014); L. Piroli, P. Calabrese, F.H.L.Essler, Phys. Rev. Lett. , 070408 (2016); L. Piroli, P.Calabrese, Phys. Rev. A , 053620 (2016). M. Rigol and A. Muramatsu Phys. Rev. Lett. , 240403(2005); M. Rigol and A. Muramatsu, Mod. Phys. Lett. B19, 861 (2005); W. Xu, M. Rigol, Phys. Rev. A , 033617(2017); A. Minguzzi and D. M. Gangardt Phys. Rev. Lett. , 240404 (2005); F. Cartarius, E. Kawasaki, and A. Min-guzzi Phys. Rev. A 92, 063605 (2015); A. del Campo andJ. G. Muga, 2006 EPL 74 M. Kormos, M. Collura, P. Calabrese, Phys. Rev. A ,013609 (2014). P. P. Mazza, M. Kormos, M. Collura, P. Calabrese, J. Stat.Mech. (2014) P11016. L. Tonks, Phys. Rev. , 955 (1936), M. Girardeau, J.Math. Phys. , 516 (1960). S. Weinberg (1996). The quantum theory of fields (Vol. 1). Cambridge university press. A. Lenard, J. Math. Phys. , , 1268-1272 (1966). E. Wigner, Phys. Rev. , 749 (1932); M. Hillery, R.F.OConnell, M.O. Scully and E.P. Wigner, Phys. Rep. (1984) 121 ; F.J. Narcowich Physica 134 A (1985) 193-208. David J. Griffiths, Introduction to quantum mechanics. Cambridge University Press (2016). K. V. Kheruntsyan, D. M. Gangardt, P. D. Drummond,and G. V. Shlyapnikov Phys. Rev. A , 053615 (2005), P.Wendenbaum, M. Collura, and D. Karevski Phys. Rev. A , 023624 (2013), Y. Brun, J. Dubail arxiv:1701.02248