Nonequilibrium quantum thermodynamics of determinantal many-body systems: Application to the Tonks-Girardeau and ideal Fermi gases
NNonequilibrium quantum thermodynamics of the Tonks-Girardeau gas
Y. Y. Atas,
1, 2
A. Safavi-Naini,
1, 3 and K. V. Kheruntsyan School of Mathematics and Physics, University of Queensland, Brisbane, Queensland 4072, Australia Institute for Quantum Computing, Department of Physics and Astronomy,University of Waterloo, Waterloo N2L 3G1, Canada ARC Centre of Excellence for Engineered Quantum Systems,School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia (Dated: May 18, 2020)We develop a general method for calculating the characteristic function of the work distribution ofthe Tonks-Girardeau (TG) gas confined to an arbitrary time-dependent trapping potential. We usethe determinant representation of the many-body wave function to characterise the nonequilibriumthermodynamics of the TG gas and obtain exact and computationally tractable expressions interms of Fredholm determinants for the mean work, the work probability distribution function, thenonadiabaticity parameter, and the Loschmidt amplitude. When applied to a harmonically trappedTG gas, our results for the mean work and the nonadiabaticity parameter reduce to those derivedpreviously using an alternative approach. We next propose to use periodic modulation of the trapfrequency in order to drive the system to highly non-equilibrium states by taking advantage of thephenomenon of parametric resonance. Under such driving protocol, the nonadiabaticity parametermay reach large values, which indicates a large amount of work being done on the system, ascompared to sudden quench protocols considered previously. This scenario is realizable in ultracoldatom experiments, aiding fundamental understanding of all thermodynamic properties of the system,in addition to facilitating possible optimization of high-efficiency quantum heat engine or quantumrefrigeration cycles.
I. INTRODUCTION
The study of thermalization in isolated quantum sys-tems has fuelled the recent development of quantumthermodynamics [1–4] and has highlighted its connectionto the dynamics of quantum correlations, scrambling ofquantum information, and entanglement dynamics [5–9].Furthermore quantum thermodynamics plays an essen-tial role in designing quantum heat engines. Here theunique quantum properties of the working fluid of the en-gine can be exploited to achieve quantum supremacy interms of the engine efficiency and output power. In thiswork we consider a system amenable to analytical andnumerical exploration in the context of quantum ther-modynamics.Recent advances in quantum simulation platforms uti-lizing a variety of atomic, molecular, and optical plat-forms has facilitated the observation of fundamental top-ics in quantum thermodynamics in the presence of differ-ent particle statistics, flexible dimensionality, and short-range and long-range interactions [10–15]. More recently,a variety of experimental platforms, including trapped-ions and nitrogen vacancy (NV) centers, have observedquantum effects in the absence of interactions [16, 17].This underscores the need for theoretical studies whichconsider the role of quantum many-body interactions.The next generation of experiments will feature quan-tum many-body interactions, allowing one to use quan-tum correlations and multi-particle entanglement as anadditional resource to improve the performance of quan-tum heat engines and quantum refrigerators. Trappedultracold atomic gases lend themselves as a particularlypromising platform in this respect [9, 10, 18–25] owing to a high degree of control over system parameters such asinter-atomic interactions and trapping potentials.The performance of a heat engine, quantum or classi-cal, is assessed using a variety of measures including theaforementioned efficiency which is given by the ratio ofthe net work done by the system and the heat energyinput during the engine cycle. The efficiency of a quan-tum heat engine requires one to calculate the mean work (cid:104) W ( t ) (cid:105) = (cid:104) ˆ H ( t ) (cid:105) − (cid:104) ˆ H (0) (cid:105) performed during the strokesof the engine cycle, in which the system is isolated fromthe heat reservoirs and hence evolves unitarily. Such evo-lution, assumed to be taking place between time zero and t in a single stroke, is described by the unitary operatorˆ U ( t ) = T e − i (cid:82) t dt (cid:48) ˆ H ( t (cid:48) ) / (cid:126) (with T being the time-orderingoperator) where the system Hamiltonian ˆ H ( t ) can in gen-eral vary in time during the evolution.The quantum nature of the system complicates the cal-culation of the mean work (cid:104) W ( t ) (cid:105) . First, W ( t ) is a ran-domly distributed quantity requiring two projective en-ergy measurements at initial time and at time t [26–30].Hence, quantum work, is not represented by a Hermi-tian operator and hence cannot be regarded as an or-dinary quantum observable [28] (see also Ref. [31, 32]).Furthermore the outcome of these projective measure-ments depend on the transition probabilities betweendifferent quantum states. This probabilistic nature im-plies that in order to fully characterize and understandthe work performed during time t one needs to knowthe full work probability distribution P ( W ) or its cor-responding characteristic function G β ( ϑ ). Finally, sincethese calculations involve many-body expectation values,their computational complexity may grow exponentiallywith system size. This is a major obstacle for a theoret- a r X i v : . [ c ond - m a t . qu a n t - g a s ] M a y ical optimization of the engine performance as it limitsany such study to small systems. In this work we showthat the Tonks-Girardeau (TG) gas [33, 34] in a time-varying potential is an ideal platform for exploring ideasof quantum thermodynamics pertaining to the design ofquantum heat machines. This system is realizable ex-perimentally [35–38], in addition to being amenable toanalytical treatment for its finite-temperature dynamics[39], and as we have shown in our previous work [40],the system dynamics in a periodically modulated har-monic trap [39, 41–43] can be stable or unstable at agiven modulation frequency depending on the amplitudeof the modulation. Hence, as we show here, we are ableto recast the quantum thermodynamics of this paradig-matic interacting many-body system in terms of finitetemperature-dynamics, which itself can be reduced to thedynamics of single-particle wavefunctions.In particular, we show that for an arbitrary modu-lation of the trapping potential thermodynamics quanti-ties of interest, such as the characteristic function G β ( ϑ ),the moments (cid:104) W n (cid:105) of the work distribution P ( W ), andthe nonadiabaticity parameter Q ∗ ( t ), can all be evalu-ated using Fredholm determinants associated solely withsingle-particle wavefunctions. Next we restrict our stud-ies to sinusoidal modulations of a harmonic trap, wherewe are able to present explicit analytical expressions forthe aforementioned quantities. We show that by sim-ply tuning the amplitude of the modulation we can sig-nificantly increase the amount of energy pumped intothe system (equivalent to the amount of work done onan isolated system) in a given time. Our expressionscan serve as a guide for designing efficient quantum en-gine cycles based on periodic modulation rather than themore traditional sudden quenches (see, e.g., [9]). Finally,our simple determinantal expression for the characteristicfunction allows us to explicitly calculate the expectationvalue (cid:104) e − βW (cid:105) (where β = 1 /k B T is the inverse tempera-ture) and hence demonstrate the validity of the quantumversion of the Jarzynski equality [27, 29, 44, 45], whichrelates quantum work associated with a nonequilibriumprocess with the free energy difference of the initial andfinal states of this quantum many-body system. II. THERMODYNAMIC QUANTITIES FORTHE TONKS-GIRARDEAU GAS IN ATIME-DEPENDENT TRAP
The Tonks-Girardeau (TG) model describes a gas of N bosons of mass m interacting in 1D via two-body hard-core interaction potential, wherein the hard-core diame-ter imposes a constraint on allowed many-body wavefunc-tions and enables the Fermi-Bose mapping [33, 34, 46–48]. Equivalently, it can be described as the limitingcase of the Lieb-Lininger model with infinitely strongtwo-body delta-function interaction potential [49]. Weassume that at time t = 0 the system is in a thermalequilibrium with an initial reservoir. It is then discon- nected from the reservoir and is subject to a confin-ing time-dependent one-body trapping potential V ( x, t )which drives the system out of equilibrium. The evolu-tion is governed by the Hamiltonianˆ H ( t ) = N (cid:88) j =1 (cid:34) − (cid:126) m ∂ ∂x j + V ( x j , t ) (cid:35) . (1)We can evaluate the work performed on the system us-ing two projective energy measurements, one performedat time t = 0 and the other at time t . The work in eachrealization of the experiment is given by the difference inthe measured energies, and the associated work proba-bility distribution function in an ensemble of realizationsis given by [27–30, 50] P ( W ) = ∞ (cid:88) N =0 (cid:88) s (cid:88) s (cid:48) p (0) N,s p ( t ) s (cid:48) | s δ ( W − E ( t ) N,s (cid:48) + E (0) N,s ) . (2)Here, p (0) N,s represents the probability of the first mea-surement returning an energy eigenvalue E (0) N,s corre-sponding to the N -particle eigenstate | Ψ s (0) (cid:105) (where weomit the index N for notational simplicity) of the ini-tial Hamiltonian ˆ H ( t = 0), satisfying ˆ H ( t = 0) | Ψ s (0) (cid:105) = E (0) N,s | Ψ s (0) (cid:105) . We consider initial thermal equilibriumstates described by the grand-canonical ensemble, andtherefore p (0) N,s is given by the normalized Gibbs factor, p (0) N,s = Z e − β ( E (0) N,s − µN ) . Here, β = 1 /k B T is the initialinverse temperature (with k B the Boltzmann constant), µ is the initial chemical potential, and Z is the corre-sponding grand-canonical partition function. The systemis then isolated from the reservoir and undergoes unitaryevolution generated by ˆ U ( t ) = T e − i (cid:82) t dt (cid:48) ˆ H ( t (cid:48) ) / (cid:126) . The re-sulting state is | Ψ s ( t ) (cid:105) = ˆ U ( t ) | Ψ s (0) (cid:105) . Next, a secondprojective energy measurement at time t returns one ofthe instantaneous energy eigenvalues E ( t ) N,s (cid:48) , with the cor-responding instantaneous eigenstate denoted via | Φ ( t ) s (cid:48) (cid:105) ,such that ˆ H ( t ) | Φ ( t ) s (cid:48) (cid:105) = E ( t ) N,s (cid:48) | Φ ( t ) s (cid:48) (cid:105) . Therefore, the sec-ond probability entering into Eq. (2), p ( t ) s (cid:48) | s , is given bythe transition probability from | Ψ s ( t ) (cid:105) to | Φ ( t ) s (cid:48) (cid:105) , namely p ( t ) s (cid:48) | s = | (cid:104) Φ ( t ) s (cid:48) | Ψ s ( t ) (cid:105) | = | (cid:104) Φ ( t ) s (cid:48) | ˆ U ( t )Ψ s (0) (cid:105) | . Finally,the delta function in Eq. (2) ensures the conservation ofenergy, such that in each realisation of the protocol thework is given by W ( t ) = E ( t ) N,s (cid:48) − E (0) N,s .However, in practice, it is often more convenient towork with the characteristic function of the work distri-bution, which is given by the Fourier transform of P ( W ), G β ( ϑ ) = (cid:90) d W e iϑW/ (cid:126) P ( W ) = (cid:104) e iϑW/ (cid:126) (cid:105) . (3)The characteristic function is the generating function ofthe moments of the distribution of work through succes-sive differentiation, (cid:104) W ( t ) n (cid:105) = ( − i (cid:126) ) n d n G β ( ϑ ) dϑ n (cid:12)(cid:12)(cid:12) ϑ =0 , (4)with the first moment corresponding to the mean work (cid:104) W ( t ) (cid:105) .Evaluating G β ( θ ) for a generic quantum many-bodysystem is a daunting task. The TG gas allows us tofind a numerically tractable expression for G β ( θ ) for ageneral V ( x, t ) using the mapping between the many-body wavefunctions Ψ s ( x , ..., x N ; t ) of hard-core bosonsand those of non-interacting fermions Ψ Fs ( x , ..., x N ; t )[33, 46–48]:Ψ s ( x , ..., x N ; t ) = A ( x , ..., x N )Ψ Fs ( x , ..., x N ; t ) . (5)with a similar relation for the instantaneous eigen-functions Φ ( t ) s ( x , ..., x N ). Here, A ( x , ..., x N ) = (cid:81) ≤ i ≤ j ≤ N sgn(x j − x i ) is the unit antisymmetric functionthat ensures the right symmetry of the bosonic wavefunc-tion under particle exchange, and Ψ Fs ( x , . . . , x N ; t ) canbe written as Slater determinants of the single-particlewavefunctions φ s i ( x, t ). Here, φ s i ( x,
0) are the energyeigenstates for the initial trapping potential V ( x,
0) witheigenenergies E (0) s i , such that the total energy is given by E (0) N,s = (cid:80) Ni =1 E (0) s i , and { s j } j =1 , ,...,N are a set of positiveinteger quantum numbers.The determinantal structure of the fermionic many-body wavefunction reduces the evaluation of Eq. (3) tothat of evaluating integrals involving time-evolved single-particle wavefunctions φ s i ( x, t ) and the instantaneoussingle-particle eigenfunctions φ ( t ) s i ( x ) with respective in-stantaneous eigenenergies E ( t ) s i . Using the properties ofFredholm integral equations and determinants [51, 52],the characteristic function can be written as (see Ap-pendix A for details), G β ( ϑ ) = D k ˆ K ( ϑ ) D f ˆ F = det(1 + ˆ K )det(1 + ˆ F ) , (6)where the superscripts k and f indicate the kernels of therespective Fredholm determinants, D k ˆ K and D f ˆ F , given by k ( x, y ) = (cid:88) i,j φ i ( x, t ) k ij ( t ) (cid:16) φ ( t ) j ( y ) (cid:17) ∗ , (7) f ( x, y ) = z (cid:88) i e − βE (0) i φ i ( x, (cid:16) φ (0) i ( y ) (cid:17) ∗ , (8)with k ij ( t ) = ze − E (0) i ( β + iϑ/ (cid:126) )+ iϑE ( t ) j / (cid:126) (cid:90) dwφ ∗ i ( ω, t ) φ ( t ) j ( ω ) , (9)and z = e βµ denoting the fugacity. Here, the operatorsˆ K and ˆ F are integral operators with kernels given by k ( x, y ) and f ( x, y ), respectively. For instance, the actionof ˆ K on an arbitrary function ξ ( r ) is given by ( ˆ Kξ )( w ) = (cid:82) R k ( w, v ) ξ ( v ) dv .Equation (6) is a significant result of this work as itallows us to express the characteristic function of thework distribution function of the quantum many-body system exclusively in terms of the single-particle quan-tities. This can be made more explicit if one rewritesEq. (6) in terms of the eigenvalues of ˆ K and ˆ F , which wedenote by { Λ ( ˆ K ) i } i =0 , ,... and { Λ ( ˆ F ) i } i =0 , ,... , respectively.The eigenvalues for the operator ˆ K (and similarly for ˆ F )and their corresponding eigenfunctions { θ ( ˆ K ) i ( w ) } i =0 , ,... satisfy the equation (cid:90) R dv k ( w, v ) θ ( ˆ K ) i ( v ) = Λ ( ˆ K ) i θ ( ˆ K ) i ( w ) , (10)with a similar equation for the operator ˆ F , where k ( w, v )is replaced by f ( w, v ). Using the expansion of the deter-minant of an operator as a product over its eigenvalues,equation (6) can be rewritten as G β ( ϑ ) = (cid:89) i (cid:32) ( ˆ K ) i ( ˆ F ) i (cid:33) , (11)which expresses the characteristic function of the workdistribution of the TG gas under an arbitrary drivingprotocol as an infinite product over the eigenvalues oftwo integral operators with kernels given by (7) and (8).While the exact analytical evaluation of the eigenval-ues is not always possible, the Fredholm determinantrepresentation offers a compact and a computationallypractical way to efficiently evaluate many-body thermo-dynamic quantities of interest numerically, in terms ofsingle-particle wavefunctions. For sufficiently smoothkernels, a simple tabulation of the kernel on a finite gridand the use of Nystr¨om classical quadrature routine en-ables one to numerically evaluate the Fredholm determi-nant with small absolute errors [52]. In Appendix B, weshow that the eigenvalues Λ ( ˆ K ) i (and Λ ( ˆ F ) i ) of the Fred-holm integral equation can also be efficiently obtainedthrough direct diagonalization of a matrix with elementsproportional to the overlap between the time evolved andinstantaneous eigenfunctions.The above discussion is written generally for a ther-mal initial state. However, in many quantum sim-ulation platforms, the relevant initial state is a purestate such as the many-body ground state of the ini-tial Hamiltonian | Ψ (0) (cid:105) . In this regime, the sensitiv-ity of the system to the driving protocol can be probedvia the Loschmidt echo, whose amplitude is equivalentto the zero-temperature characteristic function G ( t ) = G β →∞ ( t ) [53]. The Loschmidt echo amplitude is givenby (see, e.g., Ref. [54] for a review) G ( t ) = (cid:104) Ψ (0) | e i ˆ H (0) t/ (cid:126) T e − i (cid:82) t dt (cid:48) ˆ H ( t (cid:48) ) / (cid:126) | Ψ (0) (cid:105) , (12)with the Loschmidt echo itself given by L ( t ) = |G ( t ) | . (13)The Loschmidt echo gives the survival probability ofthe initial eigenstate | Ψ (0) (cid:105) first evolved forward in timeaccording to ˆ H (0) and then backward in time accord-ing to the dynamics generated by ˆ H ( t ). Its utility is incharacterising the survival of a quantum state when animperfect time-reversal is applied. Since ˆ H (0) | Ψ (0) (cid:105) = E | Ψ (0) (cid:105) and T e − i (cid:82) t dt (cid:48) ˆ H ( t (cid:48) ) / (cid:126) | Ψ (0) (cid:105) = ˆ U ( t ) | Ψ (0) (cid:105) = | Ψ ( t ) (cid:105) is the time-evolved state of the system, one canrewrite the Loschmidt echo as L ( t ) = |(cid:104) Ψ (0) | Ψ ( t ) (cid:105)| . (14)We now recognise this expression as the dynamical fi-delity, F ( t ) = |(cid:104) Ψ (0) | Ψ ( t ) (cid:105)| , which is simply thesquared overlap between the initial state and the time-evolved state of the system.In order to derive a compact and computationallytractable expression for the Loschmidt echo L ( t ), we use asimilar procedure to that used in the derivation of Eq. (6)(see Appendix E). The echo amplitude G ( t ) can then bewritten as the determinant of an N × N matrix A (where N is the number of particles in the system) containingthe overlaps between the initial and time evolved single-particle eigenfunctions φ n ( x, t ) [see Eq. (17) below], andhence L ( t ) = | det A ( t ) | , (15)where A mn ( t ) = (cid:90) ∞−∞ φ ∗ m ( x, φ n ( x, t ) dx, (16)and m, n = 0 , , . . . , N −
1. This determinantal resultreproduces the expression derived previously in Refs. [55,56].
III. TONKS-GIRARDEAU GAS IN AHARMONIC TRAP
In this section, we apply the Fredholm determinant for-malism to a time-varying harmonic potential V ( x, t ) = mω ( t ) x /
2. This model offers the advantage of beinganalytically solvable, in addition to being typical of ultra-cold atom experiments, and thus provides an ideal plat-form for studying the work distribution and other ther-modynamic quantities of interest. We point out that thethermodynamic quantities for the TG in a harmonic trapwith arbitrary time-variation of the trap frequency ω ( t )were previously derived by Jaramillo et al. in Ref. [9] us-ing the scale invariance of the corresponding many-bodywavefunction. However, our determinantal approach al-lows us to derive a more general result for the mean work (cid:104) W ( t ) (cid:105) and the nonadiabaticity parameter Q ∗ ( t ) that isvalid for an arbitrary spatial shape of V ( x, t ) beyondthe simple harmonic trapping. An immediate applica-tion of our general results to this limit—as is done inthis section—reproduces their results for the harmonictrapping.Our numerical results in the harmonic trapping regimeare further distinguished from previous studies due to the use of a nontrivial driving protocol wherein the trap po-tential is modulated sinusoidally in time with V ( x, t ) = mω ( t ) x / ω ( t ) = ω (0) [1 − α sin(Ω t )]. Here, Ω isthe modulation frequency and α characterises the mod-ulation amplitude. Under such periodic modulation, theTG trap displays a rich phase diagram in the (Ω , α ) pa-rameter space [40], including regions of stable (bounded)and unstable (exponentially growing due to parametricresonance) dynamics [42]. The latter regime is partic-ularly advantageous for creating highly nonequilibriumstates, with very large values of the nonadiabaticity pa-rameter Q ∗ ( t ) (cid:29)
1. This in turn corresponds to largeamount of work that can be done on the system accom-panied with modest changes in the trap size, which mayprove useful in optimising the performance of a quan-tum refrigerator with the TG gas serving as the workingmedium.With an arbitrary time-varying harmonic potential,the single-particle Schr¨odinger equation is exactly solv-able and the time evolved single-particle wavefunctions φ n ( x, t ) can be obtained through a simple scaling trans-formation [41, 57] given by φ n ( x, t ) = 1 √ λ φ n (cid:16) xλ , (cid:17) exp (cid:34) i mx (cid:126) ˙ λλ − i E n ( t ) t (cid:126) (cid:35) . (17)Here, the initial wavefunctions φ n ( x,
0) are the Hermite-Gauss polynomials of a quantum mechanical 1D har-monic oscillator, φ n ( x,
0) = exp (cid:0) − x / l (0) (cid:1) H n ( x/l ho (0)) π / (cid:112) n n ! l ho (0) , (18)with frequency ω (0), energy eigenvalues E (0) n = (cid:126) ω (0)( n + 1 / l ho (0) = (cid:112) (cid:126) /mω (0). Furthermore, the scaling function λ ( t ) is asolution to the Ermakov-Pinney equation [40, 58],¨ λ ( t ) + ω ( t ) λ ( t ) = ω (0) λ ( t ) , (19)with initial conditions λ (0)=1 and ˙ λ (0) = 0, whereas thetime dependent phase E n ( t ) in Eq. (17) is defined in termsof λ ( t ) and reads as E n ( t ) = (cid:126) ω (0) (cid:0) n + (cid:1) t (cid:90) t dt (cid:48) λ ( t (cid:48) ) . (20) A. The characteristic function for a thermal state
We proceed by evaluating the characteristic function G β ( ϑ ) in Eq. (11) for the general case of a thermal initialstate at inverse temperature β . The details of the deriva-tion for the eigenvalues of the Fredholm integral equa-tion involving Hermite-Gauss polynomials can be foundin Appendix C. Here we only report the final results. Fol-lowing the substitution of (B8) and (C12) in Eq. (11), thecharacteristic function of the work distribution takes theform G β ( ϑ ) = ∞ (cid:89) i =0 zξ i +1 / ze − β (cid:126) ω (0)( i +1 / , (21)where zξ i +1 / are the eigenvalues Λ ˆ Ki of the Fredholmintegral equation with kernel (A19). In Appendix C [seeEq. (C13)] we provide an explicit analytic expression for ξ which contains the dependence on ϑ . We note thatthe denominator in Eq. (21) can be recognised as theequilibrium partition function of free fermions in a har-monic trap with frequency ω (0). For ϑ = 0, we have ξ ( ϑ = 0) = e − (cid:126) βω (0) and hence G β (0) = 1, which can bedirectly seen from the definition, Eq. (3). B. The mean work
We can now use the above results to compute the meanwork performed in the system during the driving pro-tocol. Differentiating the logarithm of G β ( ϑ ) with re-spect to ϑ and using the fact that G β (0) = 1, we find G (cid:48) β ( ϑ = 0) = (log G β ( ϑ )) (cid:48) | ϑ =0 . Thus, using Eq. (21), wearrive at the expression for the mean work, (cid:104) W ( t ) (cid:105) = − i (cid:126) dG β ( ϑ ) dϑ (cid:12)(cid:12)(cid:12) ϑ =0 = (cid:18) ω ( t ) ω (0) ζ ( t ) − (cid:19) × (cid:126) ω (0) (cid:88) i z ( i + 1 / e − β (cid:126) ω (0)( i +1 / ze − β (cid:126) ω (0)( i +1 / , (22)where the time dependent parameter ζ ( t ) is given by(from Appendix C) ζ ( t ) = 12 ω (0) ω ( t ) (cid:18) ˙ λ ( t ) + ω (0) λ ( t ) + ω ( t ) λ ( t ) (cid:19) . (23)The last line in Eq. (22) corresponds to the initial ther-mal average energy of the system at inverse temperature β , i.e., (cid:104) ˆ H (0) (cid:105) = − ∂ log Z /∂β . Therefore, the meanwork performed in the system after a driving time t isgiven by (cid:104) W ( t ) (cid:105) = (cid:18) ω ( t ) ω (0) ζ ( t ) − (cid:19) (cid:104) ˆ H (0) (cid:105) , (24)and is thus proportional to the initial equilibrium internalenergy of the gas (cid:104) ˆ H (0) (cid:105) , with a (time-dependent) con-stant of proportionality given by the term in the brackets.This result provides an explicit analytic expressionthat allows one to calculate the mean work (cid:104) W ( t ) (cid:105) for awide range of Hamiltonian parameters. Furthermore, itprovides a computationally tractable expression for eval-uating the mean work done on/by the system during atime t . The only explicit unknown here is the solution to FIG. 1. The mean work (cid:104) W ( t ) (cid:105) of a driven Tonks-Girardeaugas in a harmonic trap, normalized by the initial thermal equi-librium energy (cid:104) ˆ H (0) (cid:105) . The driving protocol used here is pe-riodic modulation of the trap strength V ( x, t ) = mω ( t ) x / ω ( t ) = ω (0) [1 − α sin(Ω t )]. Panel (a) shows (cid:104) W ( t ) (cid:105) asa function of the dimensionless time Ω t and the driving fre-quency parameter a ≡ [2 ω (0) / Ω] , for the driving amplitude α = 0 .
9, whereas panel (c) is for α = 0 .
5. This parametrizationis the same as the one used in the stability diagram of Fig. 3 ofRef. [40]. Fixing the value of a (at given α ) corresponds to aspecific realisation of the driving protocol. Panels (b) and (d)show, respectively, the cuts of (cid:104) W ( t ) (cid:105) at constant a from pan-els (a) and (c), but for a longer time span. The examples of (cid:104) W ( t ) (cid:105) in (b) at a = 7 . a = 11 correspond, respectively,to stable (semi-periodic) and unstable (exponentially grow-ing) dynamics invoked by the sinusoidal drive at this valueof the modulation amplitude α = 0 .
9. Similar examples ofunstable behaviour can be realized at other values of the pa-rameter a from the high-intensity horizontal bands in panel(a), corresponding to parametric resonances at values of a slightly above a = j , where j = 1 , , , ... [40]. For compar-ison, in the examples of panel (d), which are for a smallervalue of the modulation amplitude ( α = 0 . a = 11 is no longer unstable, and as an unstable example weshow (cid:104) W ( t ) (cid:105) at the primary parametric resonance a = 1. a simple second-order ODE, the Ermakov-Pinney equa-tion (19), for the scaling parameter λ ( t ). It is worthpointing out, however, that this particular ODE affordsanalytic solutions in term of Mathieu’s functions [40], orelse can be easily solved numerically.In Fig. 1 we show the mean work (cid:104) W ( t ) (cid:105) done on/bya TG gas in a periodically modulated trap V ( x, t ) = mω ( t ) x /
2, with ω ( t ) = ω (0) [1 − α sin(Ω t )], as a func-tion of time and the dimensionless driving frequency pa-rameter a ≡ [2 ω (0) / Ω] , for two values of the modulationamplitude α . For relatively large values of α (which wenote are bounded between 0 ≤ α ≤ α = 0 .
9, we see several highintensity horizontal bands emerging dynamically withtime. These finite-width bands correspond to the valuesof a that lie within the unstable regions of the stabil-ity phase diagram of the system [40] occurring aroundand predominantly slightly above a = 1 , , , , . . . , i.e.,around integer ratios of 2 ω (0) / Ω. In these unstableregions, the energy of the system grows exponentiallywith time due to the phenomenon of parametric reso-nance. In this driving scenario, with an explicit exam-ple shown in panel (b) for a = 11, the mean work done on the system ( (cid:104) W ( t ) (cid:105) >
0) can reach very large values (cid:104) W ( t ) (cid:105) / (cid:104) ˆ H (0) (cid:105) (cid:29)
1, even though the relative change inthe instantaneous frequency ω ( t ) of the trap at time t ,i.e., at the end of any particular realisation of the driv-ing protocol, is not wildly different from ω (0); somewhatcounter-intuitively, the final value of ω ( t ) can even besmaller than ω (0), which under an adiabatic drive wouldhave to correspond to work done by the system underadiabatic expansion, corresponding to (cid:104) W ( t ) (cid:105) <
0. Incontrast, for values of a outside the parametric resonanceband, the dynamics is stable, and (cid:104) W ( t ) (cid:105) is bounded andgenerally quasiperiodic. An example of such stable dy-namics is shown in panel (b) for a = 7 .
5. We note that (cid:104) W ( t ) (cid:105) can oscillate between positive and negative val-ues.For smaller amplitudes of the drive, the widths of theparametric resonance bands along a become narrowerand only the lower order resonances get efficiently excited(see panel (c)). In panels (d) we show resonantly grow-ing (cid:104) W ( t ) (cid:105) / (cid:104) ˆ H (0) (cid:105) (cid:29) a = 1) and α = 0 .
5, which is similar to the previous ex-ample of a = 11 at α = 0 .
9. However, the dynamics for a = 11 at α = 0 . (cid:104) W ( t ) (cid:105) between negative and pos-itive values. In fact this is the behaviour that is expectedfor nearly adiabatic drive, during which the mean workalternates between being done on the system or by thesystem, following respectively the opening or tighteningthe trap as per modulation of its frequency ω ( t ). Weconfirm this below using the nonadiabaticity parameter Q ∗ ( t ), shown in Fig. 3.To summarise, the examples of Fig. 1 illustrate the richvariety of driving protocols that can be realised in a pe-riodically modulated TG gas in a harmonic trap. Thisvariety stems from the rich stability phase diagram of thesystem which contains regions of parametrically resonantunstable dynamics. Since the behaviour is fully deter-mined by two parameters, namely the driving amplitudeand the driving frequency, these examples highlight theutility of a periodically modulated TG gas as the work-ing fluid of a quantum machine. For instance it may beused to perform large irreversible work without the needfor large changes in the trap size, governed by α and themodulation frequency Ω. -3 (a)(b) FIG. 2. The work probability distribution function P ( W )versus W [in units of the initial thermal equilibrium energy (cid:104) ˆ H (0) (cid:105) ] for a TG gas in a periodically modulated harmonictrap as in Fig. 1, evaluated at time Ω t = 10. Panels (a)and (b) are, respectively, for the same parameters as the twocurves in Fig. 1 (b), representing examples of stable and un-stable dynamics. The chemical potential µ is chosen to resultin the average number of particles (cid:104) N (cid:105) = 20, and (cid:104) ˆ H (0) (cid:105) isevaluated at the temperature well below the temperature ofquantum degeneracy, k B T / (cid:104) N (cid:105) (cid:126) ω (0) = 0 . C. Work probability distribution
Given the exact analytic expression for the character-istic function of the work distribution, Eq. (21), we canevaluate not only the mean work (cid:104) W ( t ) (cid:105) of a driven TGgas, but also any higher-order moments of the work prob-ability distribution, or indeed the full probability dis-tribution P ( W ) by taking the inverse Fourier transformof Eq. (3) at a particular time instance t . In Fig. 2(a)and (b) we show representative examples of P ( W ) underthe same driving protocol as in Fig. 1, evaluated at timeΩ t = 10 for unstable and stable dynamics, respectively.In the stable regime (see Fig. 2(a)) the probabilitydistribution P ( W ) is relatively narrow and is localizedaround a relatively small values of W/ (cid:104) ˆ H (0) (cid:105) . In thiscase, the transition probabilities p s | s (cid:48) in Eq. (2) involvetransitions to instantaneous eigenstates whose eigenen-ergies are not far away from the initial eigenenergies.Therefore, the resulting transition probabilities stronglydepend on the structure of the overlaps between the par-ticular eigenstates involved. Accordingly, P ( W ) in thisexample displays a fine structure that reflects the dis-crete nature of energy levels E (0) N,s and E ( t ) N,s (cid:48) , that con-tribute to the random outcomes of the measurement of W = E ( t ) N,s (cid:48) − E (0) N,s due to sensitivity to the variousoverlaps. In contrast, in the unstable regime (see Fig.Fig. 2(b)) the distribution is broad, smooth, and centredaround relatively large values of W/ (cid:104) ˆ H (0) (cid:105) . In this case,the discrete nature of energy levels washes out as the typ-ical instantaneous energies E ( t ) N,s (cid:48) are far from the initialeigeneneriges, and all corresponding overlaps are smalland similar in magnitude.
D. The nonadiabaticity parameter
If the drive between the initial and final states is nona-diabatic in the quantum sense [59], then φ n ( x, t ) will notbe an eigenstate of the instantaneous Hamiltonian ˆ H ( t ).To quantify the degree of nonadiabaticity of the drivingprotocol, we introduce the nonadiabaticity parameter Q ∗ ( t ) = (cid:104) ˆ H ( t ) (cid:105)(cid:104) ˆ H ( t ) (cid:105) ad . (25)This parameter is the ratio of the average energy mea-sured with the actual protocol and the average energy ob-tained through an adiabatic driving. Clearly, Q ∗ ( t ) ≥ Q ∗ ( t ) = 1 when the evolution is adiabatic.For the TG gas in a harmonic trap with an arbitrarytime dependence of ω ( t ), the adiabatic driving energy isrelated to the initial energy via the relation (cid:104) ˆ H ( t ) (cid:105) ad = ω ( t ) ω (0) (cid:104) ˆ H (0) (cid:105) , (26)which can be rewritten as (cid:104) ˆ H ( t ) (cid:105) ad = (cid:104) ˆ H (0) (cid:105) /λ ( t ),where λ ad ( t ) = [ ω (0) /ω ( t )] / is the solution to theErmakov-Pinney equation (19) in the adiabatic limit¨ λ ≈
0. Equation (26) follows from the fact that (cid:104) ˆ H (0) (cid:105) = (cid:80) N,s p (0) N,s (cid:104) Ψ s (0) | ˆ H (0) | Ψ s (0) (cid:105) = (cid:80) N,s p (0) N,s E (0) N,s with E (0) N,s = (cid:126) ω (0) (cid:80) Ni =1 ( s i +1 / (cid:104) ˆ H ( t ) (cid:105) ad = (cid:80) N,s p (0) N,s E ( t ) N,s , with E ( t ) N,s = (cid:126) ω ( t ) (cid:80) Ni =1 ( s i + 1 / (cid:104) ˆ H ( t ) (cid:105) interms of (cid:104) ˆ H (0) (cid:105) , Q ∗ ( t ), ω ( t ), and ω (0). This allows us torewrite the expression for the mean work as, (cid:104) W ( t ) (cid:105) = (cid:20) ω ( t ) ω (0) Q ∗ ( t ) − (cid:21) (cid:104) ˆ H (0) (cid:105) . (27)Thus, the knowledge of (cid:104) W ( t ) (cid:105) is equivalent to theknowledge of the nonadiabaticity parameter Q ∗ ( t ) andvice versa. Indeed, by comparing Eq. (27) with Eq. (24),we identify the nonadiabaticity parameter Q ∗ ( t ) with ζ ( t ), Q ∗ ( t ) = ζ ( t ) , (28) FIG. 3. The nonadiabaticity parameter Q ∗ ( t ) for the Tonks-Girardeau gas in a periodically modulated harmonic trap. Allparameters and representative examples are the same as inFig. 1. where ζ ( t ) itself is given by Eq. (23).In Fig. 3 we show the nonadiabaticity parameter Q ∗ ( t )for a TG gas in a periodically modulated trap V ( x, t ) = mω ( t ) x /
2, with ω ( t ) = ω (0) [1 − α sin(Ω t )], as afunction of time and the driving frequency parameter a ≡ [2 ω (0) / Ω] , for the same values of α as in Fig. 1,i.e., α = 0 . α = 0 . a = 11 in panel (b) and a = 1 in panel(d) are from the unstable region. In this regime Q ∗ ( t )reaches values much greater than one, corresponding tothe driving protocol generating highly non-equilibrium fi-nal states. In contrast, at a = 11 in panel (c), the systemevolves under nearly adiabatic dynamics where Q ∗ ( t ) ≈ a = 7 . Q ∗ ( t ). E. The nonadiabaticity parameter for pure states
It is also instructive to consider driving the systemstarting from a pure state rather than a thermal state.Since the dynamics induced by the modulation of the trapfrequency is governed, in both cases, by single-particledynamics and by a single scaling parameter λ ( t ), one ex-pects to find the same Q ∗ ( t ) as the one for the thermalinitial state. For a pure initial state, the expectation val-ues appearing in Eq. (25) can be computed easily. Letsus first consider a pure state | Ψ s (cid:105) in ket notation, char-acterised by N non-negative integers s = s , s , . . . , s N .The choice of the set of integers { s i } is quite arbitrary,corresponding in general to the ground or excited many-body states, but we will see that the final result for Q ∗ ( t )is independent of this choice. The time-evolved wave-function for a pure state at time t is obtained from theSlater determinant of single-particle time evolved wave-functions,Ψ Fs ( x , . . . , x N , t ) = 1 √ N ! det ≤ n ≤ m ≤ N φ s n ( x m , t ) , (29)and has mean energy that is a simple sum of the respec-tive single-particle energy eigenvalues E s n , (cid:104) ˆ H ( t ) (cid:105) = N (cid:88) n =1 E s n ( t ) . (30)with ˆ H ( t ) φ s n ( x, t ) = E s n φ s n ( x, t ) and hence (cid:82) dx φ ∗ s n ( x, t ) ˆ H ( t ) φ s n ( x, t ) = E s n ( t ). Using Eq. (17) forthe time-evolved eigenfunctions, we find E s n ( t ) = (cid:126) ¯ ω ( t )( s n + 1 / , (31)where we have defined¯ ω ( t ) ≡ ω (0) (cid:18) ˙ λ ( t ) + ω ( t ) λ ( t ) + ω (0) λ ( t ) (cid:19) . (32)On the other hand, the adiabatic expectation valueappearing in the denominator of Eq. (25) is given by (cid:104) ˆ H ( t ) (cid:105) ad = N (cid:88) n =1 E ( t ) s n = N (cid:88) n =1 (cid:126) ω ( t ) (cid:18) s n + 12 (cid:19) . (33)Taking the ratio of these two expectation values of theHamiltonian for evaluating Q ∗ ( t ), we see that the depen-dence on the specific configuration s cancels out, and weobtain Q ∗ ( t ) = ω ( t ) ω ( t ) . (34)With ¯ ω ( t ) given by Eq. (32), we see that this expres-sion is exactly the same as the one derived previouslyfor the general case of a thermal initial state, Eq. (28).We further emphasize that the equivalence of the resultsfor Q ∗ ( t ) for a pure and thermal initial states is onlyenabled by the scale invariance of the underlying many-body dynamics, which was exploited in Ref. [9] to derivethe same expression for the nonadiabaticity parameterusing an alternative method. F. Loschmidt echo
We now turn to the discussion of the Loschmidt echoin a harmonically trapped TG gas. We remind the readerthat while all the results below are valid for an arbitrary
FIG. 4. The base l ( t ) of the Loschmidt echo L ( t ) = l ( t ) N for the Tonks-Girardeau gas in a periodically modulatedharmonic trap. All parameters and representative exam-ples are the same as in Figs. 1 and 3. The gray (dashed)line in panel (b) is the exponentially decaying prediction of l ( t ) = e πγ e − iγτ (where τ = Ω t ), with the dimensionless de-cay rate γ = | Im( ν ) | /
2, where ν is the Floquet exponent [40]equal to ν = 1 − . i in this example. For the respectivePoincare maps, explaining the typical features seen in theseLoschmidt echo curves, see Fig. 5 and text. time dependence of the trap frequency ω ( t ), the numer-ical examples will be given for the periodic modulation,consistent with the material presented earlier in this sec-tion. The exact expression of the Loschmidt echo ampli-tude, for an initial pure state with fixed number of par-ticles N in the system, can be written as (see AppendixE for details), G ( t ) = (cid:32) (cid:101) λ ( t ) λ ( t ) (cid:33) N / e i ( E (0) n −E n ( t ) ) tN / (cid:126) , (35)where E (0) n = (cid:126) ω (0) (cid:0) n + (cid:1) and with (cid:101) λ ( t ) ≡ λ ( t ) − i ˙ λ ( t ) λ ( t ) ω (0) . (36)The Loschmidt echo itself is therefore given by L ( t ) = (cid:12)(cid:12)(cid:12) (cid:101) λ ( t ) λ ( t ) (cid:12)(cid:12)(cid:12) N = l ( t ) N . (37)The Loschmidt echo, or the dynamical fidelity betweenthe initial and final states, characterises the survival orrecurrence probability of the initial state after the systemhas evolved for time t under the specific driving protocol.Under periodic modulation of the trap frequency, the dy-namics of the Loschmidt echo can be used to diagnose thestability of the TG gas for a given set of parameters. InFig. 4, we contrast the behaviour of L ( t ) in the stable andunstable regimes, for the same parameters as in Figs. 1and 3. In the stable regime, the Loschmidt echo possessesa semi-periodic behaviour since Mathieu’s function (seeRef. [40] for the mapping between the solutions to theErmakov-Pinney equation for the scaling parameter andMathieu’s function) in this region is bounded with a realand non-integer Floquet exponent ν . On the other hand,in the unstable regime the Loschmidt echo displays anexponential decay. In this region of the parameter space,the Floquet exponent ν is complex. The scaling function λ ( t ) then contains an exponentially decaying envelope,in addition to its semi-periodic behaviour [40]. Hence weobserve an overall exponential decay [see the grey dashedline in Fig. 4 (b)], with l ( t ) = e πγ e − iγτ , where the dimen-sionless decay rate is given by γ = | Im( ν ) | / τ ≡ Ω t .We can gain more insight into the behaviour of theLoschmidt echo by studying the trajectories of a classi-cal particle in a modulated harmonic trap—a strategyoften employed in semiclassical methods. The connec-tion between the TG gas in a periodically modulatedtrap and its classical counterpart can be established bynoting that the scaling function λ ( t ) which governs thedynamics of the quantum observables in our system canbe constructed by combining two independent solutions, λ ( t ) and λ ( t ) to the homogeneous version of the sameequation, i.e., Eq. (19) with the right hand side equalto zero [40]. The homogeneous equation describes themotion of a classical particle in a periodically modulatedharmonic trap, and hence the scaling function λ ( t ) is con-structed from purely classical variables or trajectories.We use Poincar´e maps to study whether the charac-teristic features observed in the behaviour of quantumobservables, such as the Loschmidt echo, manifest in thebehaviour of classical trajectories. These maps corre-spond to sampling— at regular time intervals—one of thesolutions to the classical equation of motion, say λ ( t ),in the phases space corresponding to λ ( t ) and ˙ λ ( t ).In Fig. 5 (a) and (b) we plot such Poincar´e maps rep-resenting examples of stable and unstable dynamics, re-spectively, for the same parameters as those chosen inFig. 4 (b).In the stable regime (Fig. 5 (a)) the classical trajecto-ries are confined to a finite torus. Thus, in accordance tothe Poincar´e recurrence theorem, if we allow for a suffi-ciently long evolution time, every point on the torus willbe visited at least once. This means that the startingpoint of the dynamics, λ (0) = 1 and ˙ λ (0) = 0 (or itsvicinity) is inevitably revisited after some time t , causingthe scaling solution to return to its initial value, whichcorresponds to a revival in the Loschmidt echo. We em-phasise, however, that this is an oversimplified picturebecause the general scaling solution λ ( t ) is constructedout of two independent solutions of the classical equationof motion, represented by two tori in the phase space.Consequently, an exact revival happens if a synchronousreturn to the two different respective initial points in the -1.5 -1 -0.5 0 0.5 1 1.5-2-1.5-1-0.500.511.52 -15 -10 -5 0 5 10 15-20-15-10-505101520 (a)(b) FIG. 5. The Poincare maps for stable (a) and unstable (b)dynamics. The parameter values are the same as in Fig. 4 (a), i.e. , α = 0 . a = 7 . a = 11 for panel (b). See text for for further details. phase space occurs. Furthermore, since the trajectoriesin the stable regime are confined to finite tori, this meansthat the Loschmidt echo does not decay to zero in thisregime and has a finite non zero lower-bound value.In the unstable regime (Fig. 5(b)), on the other hand,the classical trajectories spiral out (the support of thetrajectories is unbounded) and increase in magnitude,which prevents the Loschmidt echo from revivals orquasi-periodic behaviour and explains its overall decay.However, there is a transient dynamics involved here:the trajectories get trapped in a circular motion for acertain period time and then jump again into the out-ward spiral motion. The temporary circular motion cor-responds exactly to the oscillating plateaux observed inthe Loschmidt echo in the decaying example of Fig. 4. InFig. 5(b) we have highlighted the nearly circular portionof the trajectory corresponding to the plateaux region inthe Loschmidt echo in red. G. Quantum Jarzynski equality
The determinantal approach to evaluating the char-acteristic function of work distribution can also beused to explicitly verify the quantum Jarzynski equality0[27, 29, 44, 45] in the harmonically trapped TG gas withtime-varying trap frequency ω ( t ). The Jarzynski equalityrelates the change in the free energy of the system to theirreversible work along an ensemble of trajectories con-necting the initial and final states of the system. To thisend we start from the characteristic function G β ( ϑ ) for athermal initial state. Making the substitution ϑ = i (cid:126) β inEq. (C13) of Appendix C leads to ξ ( ϑ = i (cid:126) β ) = e − β (cid:126) ω ( t ) ,and therefore the expectation value (cid:104) e − βW (cid:105) , entering thequantum Jarzynski equality, can be explicitly evaluatedas (cid:104) e − βW (cid:105) ≡ G β ( ϑ = i (cid:126) β ) = ∞ (cid:89) i =0 ze − β (cid:126) ω ( t )( i +1 / ze − β (cid:126) ω (0)( i +1 / . (38)By inspecting the right hand side of Eq. (38), we seethat it is equivalent to the ratio of an instantaneous par-tition function Z of the system at an effective equilibriumtemperature β and chemical potential µ , in a trap of fre-quency ω ( t ), and the actual partition function Z of theinitial thermal equilibrium state of the system, in a trapof frequency ω (0). Therefore, Eq. (38) can be rewrittenas (cid:104) e − βW (cid:105) = ZZ ≡ e − β (Ω − Ω ) (39)where Ω = − (1 /β ) log Z and Ω = − (1 /β ) log Z denotethe grand (or Landau) potentials of the non-equilibriumand the initial equilibrium states of the system. Thegrand potentials here are given by Ω = F − µ (cid:104) N (cid:105) andΩ = F − µ (cid:104) N (cid:105) , where F and F are the respectiveHelmholtz free energies and we have used the fact thatour system evolves in the absence of heat or particle ex-change, resulting in the convservation of particle numberis in each realization of the ensemble. Hence we can sim-plify Eq. (39) to find the standard form of the aforemen-tioned Jarzynski equality given by [29, 44] (cid:104) e − βW (cid:105) = e − β ( F − F ) . (40) IV. SUMMARY
In this work we studied the dynamics of various ther-modynamic quantities of the Tonks-Girardeau gas in atime-varying potential V ( x, t ). Throughout the paperwe considered a driving protocol whereby the TG gasstarted in thermal equilibrium with a reservoir, which isreminiscent of a single stroke of a quantum heat engineor refrigerator cycle. Subsequently it was isolated fromthe reservoir and was driven out of equilibrium by V ( x, t )for a time t . Our choice of the TG gas as the workingfluid was motivated by the immediate experimental re-alizability of our proposal, as well as the analytical andnumerical tractability of the model.Our main result utilized the determinant representa-tion of the many-body wave function for deriving the characteristic function of the quantum work probabilitydistribution and all its moments in terms of the eigenval-ues of Fredholm integral operators associated exclusivelywith single-particle wavefunctions. Our results, whichmake no assumptions about the shape of the trappingpotential or the functional form of the modulation, allowfor numerically efficient evaluation of many-body ther-modynamic quantities without the need to construct thefull many-body thermal many-body state.As an immediate application of our general approach,we considered a TG gas in a harmonic trap, with anarbitrary time variation of its frequency, and presentedexplicit formulas for all the above thermodynamic quanti-ties, as well as for the nonadiabaticity parameter and theLoschmidt echo. In addition, we have validated, throughan explicit analytic derivation, the quantum Jarzynskiequality.We next elaborated on these results for a specific dy-namical protocol corresponding to a sinusoidal modula-tion of the harmonic trap strength with time. This isan experimentally relevant scenario, and our analytic re-sults lend themselves to easy exploration of the param-eter space defined by the amplitude and the frequencyof the modulation, as well as to gaining fundamentalinsights into nonequilibrium quantum thermodynamicsand its potential quantum technology applications. Fur-thermore the TG gas in a modulated harmonic trap candisplay stable or unstable dynamics at a given modula-tion frequency depending on the amplitude of the mod-ulation. Our results indicate that unstable dynamics, inwhich the system displays the phenomenon of paramet-ric resonance, has profound effect on all thermodynamicquantities. In the case of average work, this correspondsto our ability to drive the system very far from equilib-rium and hence to perform a large amount of work on thesystem in a given amount of time. Even if the system isinitialized in its zero-temperature ground state, the un-stable dynamics can be observed by the rapid decay ofthe Loschmidt echo amplitude, which characterizes thesensitivity of the system to the driving protocol and theergodicity of the semiclassical dynamics of the system.As noted above our driving protocol can be thought ofas a single stroke of the four-stroke engine cycle. Hence,our results, in particular our explicit analytic expressionsfor the harmonically trapped TG gas can serve as a guidefor designing efficient quantum engine or refrigeration cy-cles which forgo the previously studied sudden quenchprotocols in lieu of taking advantage of periodic modula-tion under the conditions of the parametric resonance. ACKNOWLEDGMENTS
This work was supported through Australian ResearchCouncil Discovery Project Grants No. DP170101423 andNo. DP190101515.1
Appendix A: Characteristic function in terms ofFredholm determinants
Here we provide the derivation of Eq. (6) of the maintext that expresses the characteristic function of the workdistribution as a ratio of two Fredholm determinants. In-serting Eq. (2) into Eq. (3), we can rewrite the charac-teristic function as G β ( ϑ ) = 1 Z (cid:88) N (cid:88) s,s (cid:48) p ( t ) s (cid:48) | s e − β ( E (0) N,s − µN )+ iϑ ( E ( t ) N,s (cid:48) − E (0) N,s ) / (cid:126) = 1 Z ∞ (cid:88) N =0 N !) (cid:88) s ,...,s N (cid:88) s (cid:48) ,...,s (cid:48) N p ( t ) s (cid:48) | s e − β ( E (0) N,s − µN )+ iϑ ( E ( t ) N,s (cid:48) − E (0) N,s ) / (cid:126) , (A1)where the conditional probability reads p ( t ) s (cid:48) | s = (cid:90) N (cid:89) i =1 dx i dx (cid:48) i Ψ s ( x , . . . , x N , t )Φ ( t ) ∗ s (cid:48) ( x , . . . , x N ) × Ψ ∗ s ( x (cid:48) , . . . , x (cid:48) N , t )Φ ( t ) s (cid:48) ( x (cid:48) , . . . , x (cid:48) N ) . (A2)In Eq. (A1), the sums over the N -particle configura-tions s and s (cid:48) have been explicitly rewritten as 2 N inde-pendent sums over the single-particle quantum numbers s , s , . . . , s N and s (cid:48) , s (cid:48) , . . . , s (cid:48) N . The many-body eigen-state energies appearing in the exponential of Eq. (A1)can be expressed in terms of the single-particle eigenen-ergies as E (0) N,s = N (cid:88) i =1 E (0) s i , (A3)and E ( t ) N,s (cid:48) = N (cid:88) i =1 E ( t ) s (cid:48) i . (A4)For a harmonically trapped TG gas, to be treated later,the single-particle energies are given explicitly by E (0) s i = (cid:126) ω (0) (cid:0) s i + (cid:1) and E ( t ) s (cid:48) i = (cid:126) ω ( t ) (cid:0) s (cid:48) i + (cid:1) , where s i =0 , , , ... and s (cid:48) i = 0 , , , ... .We now use the Fermi-Bose mapping, Eq. (5), andthe Slater determinant form of the fermionic many-bodywave functions to write the bosonic wavefunctions asΨ s ( x , . . . , x N , t ) = A ( x , ..., x N ) √ N ! det ≤ n ≤ m ≤ N φ s n ( x m , t ) , (A5)andΦ ( t ) s (cid:48) ( x , . . . , x N ) = A ( x , ..., x N ) √ N ! det ≤ n ≤ m ≤ N φ ( t ) s (cid:48) n ( x m ) , (A6)where A ( x , ..., x N ) is the unit antisymmetric functionfrom Eq. (5), whereas where φ s n ( x, t ) and φ ( t ) s (cid:48) n ( x ) are the time-evolved and the instantaneous single-particleeigenfunctions of the Hamiltonian, respectively. Insertingthese expressions back into Eq. (A2) and interchangingthe double integral with the summations over s and s (cid:48) ,we can rewrite the characteristic function as G β ( ϑ ) = 1 Z ∞ (cid:88) N =0 N !) (cid:90) N (cid:89) i =1 dx i dx (cid:48) i (cid:32)(cid:88) s Ξ s ( x , x (cid:48) ) (cid:33) × (cid:32)(cid:88) s (cid:48) ˜Ξ s (cid:48) ( x , x (cid:48) ) (cid:33) , (A7)whereΞ s ( x , x (cid:48) ) = 1 N ! (cid:32) N (cid:89) i =1 ρ s i (cid:33) det ≤ n ≤ m ≤ N φ s n ( x m , t ) × det ≤ p ≤ q ≤ N φ ∗ s p ( x (cid:48) q , t ) , (A8)and ˜Ξ s (cid:48) ( x , x (cid:48) ) = 1 N ! (cid:32) N (cid:89) i =1 ˜ ρ s (cid:48) i (cid:33) det ≤ n ≤ m ≤ N (cid:16) φ ( t ) s (cid:48) n ( x m ) (cid:17) ∗ × det ≤ p ≤ q ≤ N φ ( t ) s (cid:48) p ( x (cid:48) q ) . (A9)In Eqs. (A8) and (A9) we have used the shorthandnotation x = ( x , x , . . . , x N ), x (cid:48) = ( x (cid:48) , x (cid:48) , . . . , x (cid:48) N ) andintroduced the parameters ρ s i = √ z exp (cid:16) − E (0) s i ( β + iϑ/ (cid:126) ) (cid:17) , (A10)and ˜ ρ s i = √ z exp (cid:16) iϑE ( t ) s i / (cid:126) (cid:17) , (A11)that depend only on the instantaneous single-particle en-ergies E ( t ) s i , the inverse temperature β , and the fugacity z = e βµ of the system.The sums over functions Ξ s ( x , x (cid:48) ) and ˜Ξ s (cid:48) ( x , x (cid:48) ) ap-pearing in the integrand of Eq. (A7) have the sameform and involve products of two determinants. Let usnow show that each of these sums can be further sim-plified and written down in terms of a single determi-nant. We illustrate the derivation with (cid:80) s (cid:48) ˜Ξ s (cid:48) ( x , x (cid:48) ),and apply the Leibniz expansion to each of the two de-terminants in Eq. (A9), denoting all permutations ofthe set (1 , , . . . , N ) for each determinant via σ and ¯ σ ,respectively. Interchanging next the summation over s (cid:48) = { s (cid:48) , s (cid:48) , ..., s (cid:48) N } with the Leibniz sums over permuta-tions σ and ¯ σ , we obtain (cid:88) s (cid:48) ˜Ξ s (cid:48) ( x , x (cid:48) )= 1 N ! (cid:88) σ, ¯ σ ( − σ +¯ σ (cid:88) s (cid:48) ,...,s (cid:48) N ˜ ρ s (cid:48) · · · ˜ ρ s (cid:48) N × N (cid:89) i =1 ( φ ( t ) s (cid:48) i ( x σ i )) ∗ N (cid:89) j =1 φ ( t ) s (cid:48) j ( x (cid:48) ¯ σ j ) . (A12)2By regrouping quantities with the same indices andintroducing the instantaneous kernel g ( x, y ; t ) ≡ (cid:88) p ˜ ρ p φ ( t ) p ( x )( φ ( t ) p ( y )) ∗ , (A13)the expansion (A12) can be simplified to (cid:88) s (cid:48) ˜Ξ s (cid:48) ( x , x (cid:48) ) = 1 N ! (cid:88) σ, ¯ σ ( − σ +¯ σ N (cid:89) i =1 g ( x (cid:48) ¯ σ i , x σ i ; t ) , (A14)which can in turn be recognised as a single determinantdet ≤ n ≤ m ≤ N g ( x (cid:48) m , x n ; t ).Similarly, if we introduce the kernel f ( x, y ; t ) = (cid:88) p ρ p φ p ( x, t ) φ ∗ p ( y, t ) , (A15)then the sum (cid:80) s Ξ s ( x , x (cid:48) ) over s in (A7) is simply givenby det ≤ n ≤ m ≤ N f ( x n , x (cid:48) m ; t ).The characteristic function can therefore be written interms of a product of two determinants in the integrandof Eq. (A7): G β ( ϑ ) = 1 Z ∞ (cid:88) N =0 N !) (cid:90) N (cid:89) i =1 dx i dx (cid:48) i × det ≤ n ≤ m ≤ N f ( x n , x (cid:48) m ; t ) det ≤ n ≤ m ≤ N g ( x (cid:48) m , x n ; t ) . (A16)Thus, the original expression for the characteristic func-tion, containing a product of four determinants, has beensimplified to contain a product of two determinants, onewith the kernel (A13) and the other with the kernel(A15).We finally make use of the Andr´eief integration formula[60] (cid:90) N (cid:89) i =1 dz i det ≤ n ≤ m ≤ N A n ( z m ) × det ≤ n ≤ m ≤ N B n ( z m )= N ! det ≤ n ≤ m ≤ N (cid:18)(cid:90) dz A n ( z ) B m ( z ) (cid:19) , (A17)and eliminate the integral over the primed variable, bytaking z i = x (cid:48) i , A n ( z m ) = f ( x n , x (cid:48) m ; t ) and B n ( z m ) = g ( x (cid:48) m , x n ; t ). This gives an expression for G β ( ϑ ) that con-tains only one determinant, G β ( ϑ ) = 1 Z ∞ (cid:88) N =0 N ! (cid:90) N (cid:89) i =1 dx i det ≤ n ≤ m ≤ N k ( x n , x m ) , (A18)with the kernel k ( x, y ) = (cid:90) dwf ( x, w ; t ) g ( w, y ; t )= (cid:88) p,q φ p ( x, t ) k pq ( t, ϑ ) (cid:0) φ ( t ) q ( y ) (cid:1) ∗ , (A19) where k pq ( t, ϑ ) = ρ p ( ϑ )˜ ρ q ( ϑ ) (cid:90) dwφ ∗ p ( w, t ) φ ( t ) q ( w ) , (A20)and where we have restored the dependence on the vari-able ϑ explicitly.We note that one could have used the Andr´eief’s iden-tity and integrated over the variable x i rather than theprimed variable x (cid:48) i as the choice is arbitrary. Obviously,this alternative should not change the final result forthe characteristic function. Carrying out the integrationover x i , one would obtain the kernel k ∗ ( x, y, − ϑ ) for thecharacteristic function, which would reflect the symmetry G ∗ β ( − ϑ ) = G β ( ϑ ).Written in the form of Eq. (A18), the calculation ofthe characteristic function has so far been reduced tothe evaluation of just one determinant involving single-particle wavefunctions. However, the evaluation of themultiple integrals still constitutes a challenge such thatthis form of the characteristic function is still not practi-cal for computation.In order to further simplify our result for G β ( ϑ ) to amore tractable expression, we note that the numeratorof Eq. (A18) is amenable to a more compact and com-putational friendly form by recognising it as the minorexpansion of the Fredholm determinant belonging to thekernel k ( x, y ) [51, 52]. Finally, in order to complete theproof of our main result, Eq. (6) of the main text, itremains to show that the grand canonical partition func-tion Z ≡ Z ( t = 0) = Tr (cid:104) e − β ˆ H (cid:105) = ∞ (cid:88) N =0 (cid:88) s e − β ( E (0) N,s − µN ) , (A21)appearing in the denominator of Eq. (A18), can also bewritten as a Fredholm determinant.In order to show this, we can multiply the ex-ponential factor in Eq. (A21) by the identity 1 = (cid:82) dx dx . . . dx N | Ψ s ( x , . . . , x N , | and expand themany-body wavefunction using its Slater determinantform. Interchanging then the multiple sum over s andthe integrals, we find Z = ∞ (cid:88) N =0 N ! (cid:90) N (cid:89) i =1 dx i det ≤ n ≤ m ≤ N f ( x n , x m ) . (A22)Here, the function f ( x n , x m ) is the equilibrium kerneland corresponds to the function (A15) multiplied by √ z ,together with t = 0 and ϑ = 0; explicitly, it reads as f ( x, y ) = z ∞ (cid:88) p =0 e − iβE (0) p φ p ( x, φ ∗ p ( y, . (A23)Equation (A22) for the partition function corresponds tothe minor expansion of the Fredholm determinant be-longing to the kernel f ( x, y ).3Collecting Eqs. (A18) and (A22) together, we thus ar-rive at the final compact form of the characteristic func-tion in terms of the Fredholm determinants, Eq. (6) ofthe main text, G β ( ϑ ) = det(1 + ˆ K )det(1 + ˆ F ) , (A24)where ˆ K and ˆ F are Fredholm integral operators withkernels (A19) and (A23), respectively.Using the expansion of the two determinants inEq. (A24) in terms of products over the respective eigen-values, the characteristic function can be rewritten as G β ( ϑ ) = (cid:89) i (cid:32) ( ˆ K ) i ( ˆ F ) i (cid:33) , (A25)which is Eq. (11) of the main text. Appendix B: Numerical approach to calculating theeigenvalues of integral operators
In the general case of a TG gas in an arbitrary trappingpotential, finding the eigenvalues Λ ( ˆ K ) i of the integral op-erator with the kernel k ( x, y ) of the form (A19), which isrequired for evaluating the characteristic function (A25),has to rely on numerical approaches. (For a harmonicallytrapped TG gas, on the other hand, it possible to find an-alytic solution for this problem, and it will be presentedin Appendices C and D).In this Appendix, we outline an efficient numerical ap-proach for evaluating the eigenvalues of the Fredholmintegral equation with the kernel k ( x, y ), for situationswhen the single-particle eigenfunctions are not knownanalytically and instead have to be found numericallyas solutions to single-particle Schr¨odinger equation. Forthis, we write the Fredholm integral equation for the ker-nel k ( x, y ), Eq. (A19), (cid:90) dy k ( x, y ) θ ( ˆ K ) i ( y ) = Λ ( ˆ K ) i θ ( ˆ K ) i ( x ) , (B1)By introducing the quantities A ( i ) q = (cid:90) dy θ ( ˆ K ) i ( y ) (cid:16) φ ( t ) q ( y ) (cid:17) ∗ , (B2)we can rewrite the Fredholm integral equation as (cid:88) p,q φ p ( x, t ) k pq A ( i ) q = Λ ( ˆ K ) i θ ( ˆ K ) i ( x ) . (B3)We next multiply both side of Eq. (B3) by ( φ ( t ) m ( x )) ∗ and integrate with respect to x , yielding (cid:88) p,q c mp k pq A ( i ) q = Λ ( ˆ K ) i A ( i ) m . (B4) Here, the coefficients c mp are the coefficients of expan-sion of the time-evolved eigenstate φ p ( x, t ) in the basisof instantaneous eigenstates, φ p ( x, t ) = (cid:80) m c mp φ ( t ) m ( x ),and are given by c mp = (cid:90) dx φ p ( x, t ) (cid:16) φ ( t ) m ( x ) (cid:17) ∗ , (B5)whereas the kernel coefficients k pq are given byEq. (A20).Let now C and K be the matrices with elements( C ) ij = c ij and ( K ) ij = k ij respectively, and A ( i ) thevector with elements A ( i ) q . Then Eq. (B4) can be rewrit-ten in the matrix form,( CK ) A ( i ) = Λ ( ˆ K ) i A ( i ) . (B6)This last equation shows that the eigenvalues Λ ( ˆ K ) i ofthe Fredholm integral operator with kernel k ( x, y ) are thesame as the eigenvalues of the matrix CK , and thus it cannow be used for numerical implementation of our eigen-value problem. In practice, one first calculates the matrix C of the overlaps between the time-evolved and instan-taneous eigenstates. The elements of the matrix K arethen easily obtained through the relation ( K ) ij = ρ i ˜ ρ j c ∗ ji ,where ρ i and ˜ ρ j are given by Eqs (A10) and (A11). Thematrix CK is finally constructed and diagonalized with asufficiently large size to ensure the convergence of eigen-values.For kernels of the form f ( x, y ) and g ( x, y ), Eqs. (A15)and (A13), the eigenvalue problem is much simpler. Suchkernels are known as degenerate kernels in Fredholm the-ory, and the eigenvalues of the associated Fredholm inte-gral equation are, in fact, given by ρ i and ˜ ρ i , respectively,without the need for any further calculations. In orderto prove this, one can follow the same argument as theone above for the kernel k ( x, y ) and show that the ma-trix equation obtained is already diagonal. This specialproperty is due to the orthogonality of the eigenfunctionsappearing in the expansions (A13) and (A15).In particular, for the degenerate kernel f ( x, y ),Eq. (A23), which appears in the determinantal form ofthe partition function Z , Eq. (A22), we observe thatit can be expressed in terms of the kernel f ( x, y ) as f ( x, y ) ≡ √ zf ( x, y, t = 0 , ϑ = 0). Therefore, its eigen-values Λ ( ˆ F ) k , required for the evaluation of the denomi-nator of the characteristic function (A25), are easily ob-tained (with the use of Eq. (A10)) asΛ ( ˆ F ) k = √ zρ k ( ϑ = 0) = z exp( − βE (0) k ) . (B7)With this results, we thus arrive at the familiar form ofthe grand-canonical partition function, Z = ∞ (cid:89) i =0 (1 + Λ ( ˆ F ) i ) = ∞ (cid:89) i =0 (1 + ze − βE (0) i ) , (B8)describing free fermions in 1D, which is the same as theone for the TG gas due to Bose-Fermi mapping.4 Appendix C: Derivation of the characteristicfunction of a harmonically trapped TG gas for athermal state
In this appendix, we show that for a harmonicallytrapped TG gas the eigenvalues Λ ( ˆ K ) i of the integral op-erator with the kernel k ( x, y ), appearing in the charac-teristic function (A25), can be calculated analytically.Indeed, for a harmonically trapped TG gas, we make useof the Hermite-Gauss polynomials for the time-evolvedand instantaneous single-particle wavefunctions in theexpressions for the kernels g ( x, y ) and f ( x, y ), given byEqs. (A13) and (A15). The kernel k ( x, y ) can then befound from g ( x, y ) and f ( x, y ), according to Eq. (A19).In order to show this, we first observe that the kernels g ( x, y ) and f ( x, y ) have the following generic functionalform, ae − b ( x + y ) ∞ (cid:88) p =0 H p ( cx ) H p ( cy )2 p p ! d p , (C1)where the constants a, b, c and d (independent of x , butdependent on time t ) are known and depend on the pa-rameters of the single-particle wavefunctions. Owing toMehler’s summation formula [61], ∞ (cid:88) p =0 H p ( cx ) H p ( cy )2 p p ! d p =(1 − d ) − / exp (cid:18) xyd − ( x + y ) d c (1 − d ) (cid:19) , (C2)the kernels g ( x, y ) and f ( x, y ) can be simplified to ageneric form of a Gaussian quadratic in x and y withknown coefficients.The Gaussian quadratic form for the the kernels g ( x, y )and f ( x, y ) allows one to also simplify the expression forthe kernel k ( x, y ), Eq. (A19), which depends on the prod-uct of g ( x, y ) and f ( x, y ). Specifically, one can integratethe product of the two Gaussian quadratic forms to ob-tain another Gaussian quadratic, k ( x, y ) = A exp (cid:0) − Bx − Cy + Dxy (cid:1) , (C3)where the coefficients A, B, C and D (which are time de- pendent) are given by A = zλl ho (0) (cid:115) uvπκ (1 − u )(1 − v ) , (C4) B = 12 λ l (0) (cid:32) u − u − u λ ad κ (1 − u ) λ − i ˙ λλω (0) (cid:33) , (C5) C = 12 l (0) λ ad (cid:18) v − v − v κ (1 − v ) (cid:19) , (C6) D = 4 uvl (0) λ κ (1 − u )(1 − v ) . (C7)Here, l ho (0) = ( (cid:126) /mω (0)) / is the harmonic oscillatorlength for the initial trap frequency ω (0), λ ( t ) is thescaling solution, λ ad ( t ) = ( ω (0) /ω ( t )) / is the scalingsolution in the adiabatic limit, z = e βµ is the fugacity,and the quantities u , v , and κ are defined according to u = e − (cid:126) ω (0)( β + iϑ/ (cid:126) ) , (C8) v = e iϑω ( t ) , (C9) κ = (cid:18) v − v (cid:19) + (cid:18) λ ad λ (cid:19) (cid:32) u − u + i ˙ λλω (0) (cid:33) . (C10)We momentarily pause here in order to refer the readerto Appendix D, in which we derive the eigenvalues Λ ( ˆ K ) i of the integral operator with the kernel k ( x, y ) given inthe general form of a Gaussian quadratic (C3). The finalresult is expressed only in terms of the coefficients A, B, C and D of the quadratic, and reads as [from Eq. (D11)]Λ ( ˆ K ) i = √ πAD i (cid:104) B + C + (cid:112) ( B + C ) − D (cid:105) i +1 / . (C11)Substituting now the expressions for the coefficients A, B, C and D from Eqs. (C4)–(C7) into Eq. (C11), aftera little algebra, we finally obtainΛ ( ˆ K ) i = zξ i +1 / , (C12)where we have introduced ξ = 4 uv (1 + u )(1 + v ) + (1 − u )(1 − v ) ζ ( t ) + (cid:113) [(1 + u )(1 + v ) + (1 − u )(1 − v ) ζ ( t )] − u v , (C13)and where ζ ( t ) = 12 ω (0) ω ( t ) (cid:18) ˙ λ ( t ) + ω (0) λ ( t ) + ω ( t ) λ ( t ) (cid:19) . (C14)In the main text, we show that the nonadiabaticity pa- rameter Q ∗ ( t ) is equivalent to the above ζ ( t ).Combining Eq. (C12) with the expression (B8) for thegrand-canonical partition function, we can now rewriteEq. (11) of the main text (or equivalently Eq. (A25)) for5the characteristic function as G β ( ϑ ) = ∞ (cid:89) i =0 zξ i +1 / ze − βE (0) i , (C15)which is Eq. (21) of the main text. Appendix D: Eigenspectrum of a Gaussian quadratickernel
As shown in the previous appendix, for a harmonicallytrapped TG gas, the different kernels (A13), (A15) and(A19) reduce to a Gaussian quadratic form. Asn an ex-ample, the kernel k ( x, y ) can be parametrized as k ( x, y ) = A exp (cid:0) − Bx − Cy + Dxy (cid:1) , (D1)where A, B, C and D are known coefficients.In this appendix, we derive an exact expression forthe eigenfunctions and eigenvalues of a Fredholm integralequation (cid:90) R dv k ( x, y ) θ ( ˆ K ) i ( y ) = Λ ( ˆ K ) i θ ( ˆ K ) i ( x ) (D2)with the Gaussian kernel k ( x, y ).We seek the solution for the eigenfunctions θ ( ˆ K ) i ( x ) inthe form of θ ( ˆ K ) i ( x ) = γ exp (cid:0) − ηx (cid:1) H i ( ζx ) , (D3)where H i ( x ) is the Hermite polynomial of order i , γ isfixed by the normalisation, whereas η and ζ are constantschosen to satisfy the integral equation (D2).Inserting the ansatz (D3) in the integral equation (D2),and making use of the formula 7 . (cid:90) R e − ( x − y ) H i ( αx ) dx = √ π (1 − α ) i/ H i (cid:18) αy (1 − α ) / (cid:19) , (D4)we find that the left hand side of Eq. (D2) is given by Aγ √ π √ C + η exp (cid:18) − w (cid:18) B − D C + η ) (cid:19)(cid:19) (cid:18) − ζ ( C + η ) (cid:19) i/ × H i (cid:32) Dζw (cid:112) ( C + η ) − ζ ( C + η ) (cid:33) . (D5)Therefore, by identification with the right hand side of(D2), we deduce the following identities: η = B − D C + η ) , (D6)1 = D (cid:112) ( C + η ) − ζ ( C + η ) , (D7)Λ ( ˆ K ) i = A √ π √ C + η (cid:18) − ζ ( C + η ) (cid:19) i/ . (D8) Equation (D8) defines the eigenvalues of the problemin terms of ζ and η . Solving the quadratic equation (D6)for η gives η = 12 (cid:16) B − C ± (cid:112) ( B + C ) − D (cid:17) , (D9)where one must choose the positive branch such thatRe( η ) > ζ = ( C + η ) (cid:18) − D C + η ) (cid:19) , (D10)which leads to the following final form for the eigenvalues,Λ ( ˆ K ) i = A √ πD i i ( C + η ) i +1 / = √ πAD i (cid:104) B + C + (cid:112) ( B + C ) − D (cid:105) i +1 / , (D11)expressed only in terms of the parameters A, B, C and D of the original Gaussian kernel (D1). We note that theeigenvalues are invariant under the exchange B ↔ C ,which means that the integration in Eq. (D2) can bedone over either the first or second variable of the kernel k ( x, y ). This in turn means that the kernels k ( x, y ) and k ( y, x ) have the same eigenvalue spectrum. Appendix E: Determinantal form of the Loschmidtecho
The Loschmidt echo (14) is given by the squared over-lap between the time-evolved many-body ground state | Ψ ( t ) (cid:105) , and the initial one at time t = 0, | Ψ (0) (cid:105) . For N particles in the ground state at t = 0, the time-evolvedand initial states are constructed with the Slater determi-nant of the corresponding N first single-particle orbitalsat time t and time t = 0, respectively. Explicitly, we haveΨ ( x , ..., x N ; t ) = A ( x , ..., x N ) √ N ! det ≤ n ≤ m ≤ N φ n − ( x m , t ) , (E1)where A ( x , ..., x N ) is the unit antisymmetric functionfrom Eq. (5). Similarly, setting t = 0 in Eq. (E1), givesthe initial ground-state wavefunction. The overlap be-tween the two wavefunctions is thus obtained as an N -fold integral over the product of two determinants (cid:104) Ψ (0) | Ψ ( t ) (cid:105) = 1 N ! (cid:90) N (cid:89) i =1 dx i det ≤ n ≤ m ≤ N φ ∗ n − ( x m , × det ≤ n ≤ m ≤ N φ n − ( x m , t ) . (E2)Owing to Andreief’s integration formula (A17), thisexpression can be further reduced to a form containing6just one determinant, (cid:104) Ψ (0) | Ψ ( t ) (cid:105) = det ≤ n ≤ m ≤ N − (cid:90) dx φ ∗ m ( x, φ n ( x, t ) . (E3)Introducing a N × N matrix A , with the matrix ele-ments corresponding to the overlaps between the time-evolved and initial single-particle wavefunctions, A mn ( t ) = (cid:90) dx φ ∗ m ( x, φ n ( x, t ) (E4)with m, n = 0 , , . . . , N −
1, one obtains the Loschmidtecho in a simple determinantal form: L ( t ) = | det A | . (E5)Evaluating the matrix elements of A ( t ) for a harmon-ically trapped TG gas, with Hermite-Gauss polynomialsfor the single-particle eigenfunctions, Eq. (18), and thetime-evolved wavefunctions, Eq. (20), gives A mn ( t ) = (cid:0) m + n m ! n ! πλ ( t ) (cid:1) − / e − i E n ( t ) t/ (cid:126) J mn ( t ) . (E6)Here, E n ( t ) is given by Eq. (20) of the main text, whereasthe matrix elements J mn ( t ) are given by J mn ( t ) = (cid:90) ∞−∞ H m ( y ) H n (cid:18) yλ ( t ) (cid:19) e − (cid:101) λ ( t ) y / dy, (E7) where y = x/l ho (0) is the dimensionless coordinate and˜ λ ( t ) is given by Eq. (36) of the main text.Since the parity of the wave function is conserved dur-ing the time evolution, transitions between single-particleeigenfunctions with different parity are not allowed. Letus define the k -th diagonal of a matrix such that k = 0corresponds to the diagonal, k = ± J mn consequently has a special alternating bandstructure with the k -th diagonal equal to zero for odd k .Using the generating function method, we find that thematrix elements J mn ( t ) are given by J mn ( t ) = (cid:114) π (cid:101) λ ¯ n ! i ¯ n ( m + n ) / λ n (cid:32) (cid:101) λ − (cid:101) λλ − (cid:33) ( m − n ) / × ( (cid:101) λλ − b − (cid:101) λ ( m + n ) / P | m − n | / m + n ) / (cid:16) λ | (cid:101) λ | (cid:17) , (E8)where ¯ n = min( m, n ) and P ql ( z ) are the associated Leg-endre polynomials.At first sight, it does not seem obvious that the ma-trix A ( t ) with such matrix elements will have a simpledeterminant as to result in the Loschmidt echo ampli-tude G ( t ) = det A ( t ) given by Eq (35) and hence thefinal results for the Loschmidt echo itself, Eq. (37), for N particles in the system. However, it is easy to guessthe exact general formula by looking at the determinantof this matrix (of size N × N ) for the first few values of N and one finds the formula (37) of the main text. Ourresult for G ( t ) agrees with that of Ref. [63] derived usingan alternative approach. [1] J. Gemmer, M. Michel, and G. Mahler, Quantumthermodynamics: Emergence of thermodynamic behaviorwithin composite quantum systems , Vol. 784 (Springer,2009).[2] F. Binder, L. A. Correa, C. Gogolin, J. Anders, andG. Adesso,
Thermodynamics in the quantum regime (Springer, 2018).[3] S. Vinjanampathy and J. Anders, Contemporary Physics , 545 (2016).[4] R. Kosloff, Entropy , 21002128 (2013).[5] S. W. Kim, T. Sagawa, S. De Liberato, and M. Ueda,Phys. Rev. Lett. , 070401 (2011).[6] J. M. Diaz de la Cruz and M. A. Martin-Delgado, Phys.Rev. A , 032327 (2014).[7] R. Nandkishore and D. A. Huse, Annu. Rev. Condens.Matter Phys. , 15 (2015).[8] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol,Adv. Phys. , 239 (2016).[9] J. Jaramillo, M. Beau, and A. del Campo, New Journalof Physics , 075019 (2016).[10] Y. Zheng and D. Poletti, Phys. Rev. E , 012110 (2015).[11] M. G¨arttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall,J. J. Bollinger, and A. M. Rey, Nature Physics , 781(2017). [12] T. Schweigler, V. Kasper, S. Erne, I. Mazets, B. Rauer,F. Cataldini, T. Langen, T. Gasenzer, J. Berges, andJ. Schmiedmayer, Nature , 323 (2017).[13] N. Friis, O. Marty, C. Maier, C. Hempel, M. Holz¨apfel,P. Jurcevic, M. B. Plenio, M. Huber, C. Roos, R. Blatt,and B. Lanyon, Phys. Rev. X , 021012 (2018).[14] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch,C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F.Roos, Science , 260 (2019).[15] K. A. Landsman, C. Figgatt, T. Schuster, N. M. Linke,B. Yoshida, N. Y. Yao, and C. Monroe, Nature , 61(2019).[16] D. von Lindenfels, O. Gr¨ab, C. T. Schmiegelow,V. Kaushal, J. Schulz, M. T. Mitchison, J. Goold,F. Schmidt-Kaler, and U. G. Poschinger, Phys. Rev.Lett. , 080602 (2019).[17] J. Klatzow, J. N. Becker, P. M. Ledingham, C. Weinzetl,K. T. Kaczmarek, D. J. Saunders, J. Nunn, I. A. Walm-sley, R. Uzdin, and E. Poem, Phys. Rev. Lett. ,110601 (2019).[18] S. Deffner and E. Lutz, Phys. Rev. E , 021128 (2008).[19] O. Abah and E. Lutz, EPL (Europhysics Letters) ,60002 (2016).[20] S. Sotiriadis, A. Gambassi, and A. Silva, Phys. Rev. E , 052129 (2013).[21] G. D. Chiara, A. J. Roncaglia, and J. P. Paz, New Jour-nal of Physics , 035004 (2015).[22] M. Beau, J. Jaramillo, and A. Del Campo, Entropy ,168 (2016).[23] C. Rylands and N. Andrei, Phys. Rev. B , 064308(2019).[24] G. Perfetto, L. Piroli, and A. Gambassi, Phys. Rev. E , 032114 (2019).[25] J. Li, T. Fogarty, S. Campbell, X. Chen, and T. Busch,New Journal of Physics , 015005 (2018).[26] S. Vinjanampathy and J. Anders, Contemporary Physics , 545 (2016).[27] J. Kurchan, arXiv:cond-mat/0007360 (2000).[28] P. Talkner, E. Lutz, and P. H¨anggi, Phys. Rev. E ,050102 (2007).[29] J. Yi, Y. W. Kim, and P. Talkner, Phys. Rev. E ,051107 (2012).[30] M. Campisi, P. H¨anggi, and P. Talkner, Rev. Mod. Phys. , 771 (2011).[31] A. J. Roncaglia, F. Cerisola, and J. P. Paz, Phys. Rev.Lett. , 250601 (2014).[32] F. Cerisola, Y. Margalit, S. Machluf, A. J. Roncaglia,J. P. Paz, and R. Folman, Nat. Commun. , 1 (2017).[33] M. Girardeau, Journal of Mathematical Physics , 516(1960).[34] M. D. Girardeau, Phys. Rev. , B500 (1965).[35] B. Paredes, A. Widera, V. Murg, O. Mandel, S. F¨olling,I. Cirac, G. V. Shlyapnikov, T. W. H¨ansch, and I. Bloch,Nature , 277 (2004).[36] T. Kinoshita, T. Wenger, and D. S. Weiss, Science ,1125 (2004).[37] E. Haller, M. Gustavsson, M. J. Mark, J. G. Danzl,R. Hart, G. Pupillo, and H.-C. N¨agerl, Science ,1224 (2009).[38] J. M. Wilson, N. Malvania, Y. Le, Y. Zhang, M. Rigol,and D. S. Weiss, Science , 1461 (2020).[39] Y. Y. Atas, D. M. Gangardt, I. Bouchoule, and K. V.Kheruntsyan, Phys. Rev. A , 043622 (2017).[40] Y. Y. Atas, S. A. Simmons, and K. V. Kheruntsyan,Phys. Rev. A , 043602 (2019).[41] A. Minguzzi and D. M. Gangardt, Phys. Rev. Lett. ,240404 (2005).[42] E. Quinn and M. Haque, Phys. Rev. A , 053609 (2014).[43] Y. Y. Atas, I. Bouchoule, D. M. Gangardt, and K. V.Kheruntsyan, Phys. Rev. A , 041605(R) (2017). [44] C. Jarzynski, Phys. Rev. Lett. , 2690 (1997).[45] M. Esposito, U. Harbola, and S. Mukamel, Rev. Mod.Phys. , 1665 (2009).[46] M. D. Girardeau and E. M. Wright, Phys. Rev. Lett. ,5691 (2000).[47] R. Pezer and H. Buljan, Phys. Rev. Lett. , 240403(2007).[48] V. Yukalov and M. Girardeau, Laser Physics Letters ,375 (2005).[49] E. H. Lieb and W. Liniger, Phys. Rev. , 1605 (1963).[50] C. Jarzynski, H. T. Quan, and S. Rahav, Phys. Rev. X , 031038 (2015).[51] A. Lenard, Journal of Mathematical Physics , 1268(1966).[52] F. Bornemann, Mathematics of Computation , 871(2010).[53] A. Silva, Phys. Rev. Lett. , 120603 (2008).[54] T. Gorin, T. Prosen, T. H. Seligman, and M. ˇZnidariˇc,Physics Reports , 33 (2006).[55] K. Lelas, T. ˇSeva, and H. Buljan, Phys. Rev. A ,063601 (2011).[56] K. Lelas, T. ˇSeva, H. Buljan, and J. Goold, Phys. Rev.A , 033620 (2012).[57] A. Perelomov and Y. Zel’dovich, Quantum Mechanics (World Scientific, Singapore, 1998).[58] V. P. Ermakov, Univ. Izv. Kiev , 1 (1880).[59] We note here that the term adiabatic in a thermody-namic sense is often used to merely imply thermal iso-lation from the environment and hence absence of anyheat transfer between the system and the environment.Within this definition, an adiabatic process can be re-versible (in which case the process would be isentropic, i.e. , entropy preserving) or irreversible (entropy generat-ing). In contrast, adiabatic in the quantum sense refers tosatisfying the quantum adiabatic theorem, and thereforea nonadiabatic evolution of an otherwise thermally iso-lated system implies that, during the evolution, quantumtransitions may occur in general, and hence the processis irreversible and is entropy generating.[60] P. Andr´eief, Mm. Soc. Sci. Phys. Nat. Bordeaux , 1(1886).[61] H. Bateman, Higher Transcendental Functions [VolumesI-III] , Vol. 3 (McGraw-Hill Book Company, 1953).[62] I. S. Gradshteyn and I. M. Ryzhik,
Table of integrals,series, and products (7th. Ed., Academic Rress, 2017).[63] A. del Campo, New Journal of Physics18