Functional field integral approach to the quantum work

# Functional field integral approach to the quantum work

Jian-Jun Dong Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences, Beijing 100049, China    Yi-feng Yang Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences, Beijing 100049, China Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China Collaborative Innovation Center of Quantum Matter, Beijing 100190, China
July 20, 2019
###### Abstract

We introduce the functional field integral approach to study the statistics of quantum work for arbitrary nonequilibrium process. For a general bilinear Hamiltonian with a sudden quench of the control parameters, we derive the analytical formalism for the characteristic function, the mean work density and mean irreversible work density. We show that the mean irreversible work may be used to detect quantum phase transitions and exhibit characteristic scaling at the quantum critical point for a single quench process. For a double quench process, a dynamical quantum phase transition may occur as a singularity in the so-called rate function. We examine our method in three models: the transverse Ising chain, the Su-Schrieffer-Heeger (SSH) model and the Bardeen-Cooper-Schrieffer (BCS) Hamiltonian for superconductivity. Our results for the Ising chain is in agreement with the Hamiltonian approach and yield the correct scaling exponents, while for the SSH model, we find nonuniversal quantum critical scaling due to the topological nature of its quantum phase transition. Quantum critical scaling and dynamical quantum phase transition were also obtained in the BCS model with single or double quench protocols in combination with the mean-field approximation. This suggests that our approach may be potentially extended to other many-body systems for future investigation of their dynamical properties under nonequilibrium condition.

## I Introduction

Nonequilibrium conditions provide an additional dimension in time domain for probing the many-body dynamics beyond the well-established equilibrium statistics and have recently led to the proposal of many novel phenomena such as the dynamical quantum phase transition Heyl2013 (); Budich2016 (); Heyl2018 (); Lang2018 (), the time crystal Wilczek2012 (); Else2016 (), the fluctuation relations Talkner2007PRE (); Talkner2007JPA (); Crooks1999 (); Tasaki2000 (); Jarzynski1997PRL (); Jarzynski1997PRE (); Esposito2009 (); Campisi2011 (), and so on. While ultracold atoms in optical lattices can be easily tuned to be out of equilibrium Bloch2008 (); Aoki2014 (), recent development of ultrafast pump-probe spectroscopy has enabled one to study the excitation and relaxation dynamics in real correlated materials Wall2011 (). The study of nonequilibrium quantum physics is becoming one of the most active and exciting branches of modern condensed matter physics and attracted intensive attentions in recent years Polkovnikov2011 (). However, despite of many theoretical progresses including the nonequilibrium extension of the density matrix renormalization group, quantum master equations, Keldysh Green’s function technique and dynamical mean-field theory Schollwock2005 (); Kamenev2009 (); Breuer2002 (), we still lack a schematic framework to interpret the vast kinds of nonequilibrium phenomena.

Probability distribution of quantum work is arguably one of the most important quantities to characterize the nonequilibrium dynamics such as the dynamical fluctuations, quantum phase transitions (QPTs) and quantum criticality Imparato2005 (); Engel2007 (); Palmai2014 (); Fusco2014 (); Talkner2016 (); Lobejko2017 (); Modak2017 (); Wang2017 (); Funo2017 (); Russomanno2015 (); Shraddha2015 (); Jarzynski2015 (); Nigro2018 (). Technically, it may be defined as

 p(w)=∑n,mδ(W−Efm+Ein)P(mf∣ni)P(ni), (1)

where is the work density, is the lattice size, is the dimensionality, is the -th eigenvalue of the initial state (), and is the -th eigenvalue of the final state (). denotes the probability distribution of the initial canonical state, and accounts for the transition probability between the initial and final states governed by the time-dependent evolution, , where is the time ordering operator and is the period for the time evolution. The reduced Planck constant is set to unity. The nonequilibrium dynamics is fully incorporated in the time-dependent Hamiltonian and the transition probability contains the key information characterizing the dynamical process.

The above method has been applied to a variety of many-body systems including the transverse Ising chain Silva2008 (); Dorner2012 (), the anisotropic XY model Bayocboc2015 (), the XXZ model Mascarenhas2014 (), the Luttinger liquid Dora2012 (), and the low-dimensional quantum gas Gambassi2012 (); Sotiriadis2013 (); Shchadilova2014 (). In most of these works, a Hamiltonian approach has been used by calculating explicitly the time evolution of the quantum state for quench protocol. However, such calculations can be very involved which makes it impossible to be extended to more general correlated many-body systems Dora2012 (); Smacchia2013 (). To overcome this issue, here we propose to carry out the calculations straightforwardly using the functional field integral formalism. The latter has the advantage of numerical simplicity and may be extended to more complicated nonequilibrium processes beyond the quench protocol. Moreover, in combination with the auxiliary field method and mean-field approximation, it has the potential to be applied to more general correlated systems, in the cost of explicit state evolution.

The manuscript is organized as follows. In Section II, we first introduce the general formalism based on the functional field integral approach and show that it can be used to identify both the equilibrium and dynamical QPTs. In Section III, we apply it to three different models, including the well-studied transverse Ising model, the Su-Schrieffer-Heeger (SSH) model with topological phase transition and the Bardeen-Cooper-Schrieffer (BCS) model with superconductivity under mean-field approximation. We will discuss the QPTs in these models and derive the corresponding quantum critical exponents from the calculated mean irreversible work density. Section IV is the discussion and conclusions.

## Ii General formalism

To proceed, we consider the bilinear model,

 H(t)=∑kΨ†kAk(t)Ψk, (2)

where is a matrix and is the vector of the Pauli matrices. Although simple, it represents a large number of models in condensed matter physics. We first reformulate the functional field integral approach to calculate the Fourier transformation of the work distribution, namely the characteristic function Talkner2007PRE ():

 G(u)=∫∞−∞dwexp(iuwNd)p(w) (3)

where is the partition function at . The time contour of the functional field integral in is illustrated in Fig. 1. After some tedious calculations using the functional field integral techniques (see Appendix A), we obtain

 (4)

where

 Bk(T0)=C†k(T0)eiuAk(T0)Ck(T0)e−(iu+β)Ak(0), (5)

with . The mean work density is given by the first cumulant of the characteristic function, , yielding

 ⟨w⟩=1Nd∑kTr[(Dk(T0)−Ak(0))e−βAk(0)]2+Tr[e−βAk(0)], (6)

with , where can be computed numerically for arbitrary time dependence. Note that the mean work density always exceeds the free energy density difference between the initial and final equilibrium states (both with the same inverse temperature ), as stated in the second law of thermodynamics. A mean irreversible work density can thus be defined as their difference, , which characterizes the irreversibility of the nonequilibrium process and is directly related to the entropy increase between the final and initial equilibrium states, , for a closed quantum system without heat transfer Shraddha2015 ().

The mean work density and the mean irreversible work density can be used to identify equilibrium phase transitions. Considering a quench process where the model Hamiltonian changes from at time to at , we have . As shown in Appendix B, one can derive an explicit formula for the free energy density difference,

 Δf=−1βNd∑klncosh2(βE1k/2)cosh2(βE0k/2), (7)

and the mean work density,

 ⟨w⟩=1Nd∑k[(E0k)2−d0k⋅d1k]tanh(12βE0k)E0k, (8)

where . Then the mean irreversible work density can be immediately obtained by definition. To see how these detect the phase transitions, we consider a small quench, . It can be shown analytically (Appendix B) that and at zero temperature, which connect directly to the first and second derivatives of the ground state energy with respect to the quench parameter. Thus the singularity in and reflect the first or second-order phase transitions.

For a (second-order) QPT, scaling analysis of can yield key information on the critical exponents Shraddha2015 (). In the heat susceptibility limit, where is the largest length scale, the mean irreversible work density at zero temperature scales as

 ⟨wirr⟩/δ2 ∼λν(d+z)−2,  δ−ν>N>λ−ν ∼N2/ν−(d+z),δ−ν>λ−ν>N, (9)

where is the correlation length exponent, is the dynamical exponent, and reflects the distance from the quantum critical point (QCP). In the thermodynamic limit, where is the largest length scale, one arrives at the scaling relation

 ⟨wirr⟩/δ2 ∼δν(d+z)−2,N>λ−ν>δ−ν ∼λν(d+z)−2,N>δ−ν>λ−ν. (10)

On the other hand, if the system is prepared very close to the QCP such that is larger than all other length scales, the scaling relation becomes

 ⟨wirr⟩/δ2 ∼δν(d+z)−2,  λ−ν>N>δ−ν ∼N2/ν−(d+z),  λ−ν>δ−ν>N. (11)

Note that when the combination exceeds 2, the scaling of is non-universal. In the marginal case, , one may find additional logarithmic correction to above scalings Shraddha2015 ().

Recently, it has also been shown that the work statistics in a double quench process, where the Hamiltonian changes from to at and then back to for , can be used to describe dynamical QPTs Heyl2013 (). In this case, the free energy density difference is zero and we have and . The mean irreversible work density is incapable of capturing the dynamical QPTs. On the other hand, the work distribution function, which in principle could be evaluated by inverse Fourier transformation of , is often ill-defined in numerical calculations. Lately, an alternative approach has been proposed based on the so-called Gärter-Ellis theorem Abeling2016 (); Majumdar2017 (). Considering a global quench process, the quantum work grows exponentially with the system size, , which defines the rate function . The theorem states that can be obtained via a Legendre-Fenchel transformation,

 r(w)=−infR∈R[wR−c(R)], (12)

where is a scaled cumulant generating function assumed to be differentiable with respect to the real variable , and the infimum is evaluated within the domain of definition of including . The dynamical QPT is then manifested as a singular point of , as discussion in Appendix B.

## Iii Numerical results

In this section, we use the above general formalism to study the work statistics and possible QPTs in three different models. The first model is the transverse Ising model, which has been extensively studied using the Hamiltonian approach and hence provides a good examination of our method. We then discuss the SSH model, where we will find corrections to the quantum work due to its topological properties. Last but not least, for possible extension to other correlated models in future studies, we discuss the well-known BCS model for superconductivity in combination with the mean-field approximation.

### iii.1 The transverse Ising model

The Hamiltonian of the transverse Ising model is Pfeuty1970 (); Suzuki2013 (); Dutta2015 ()

 H=−JN∑j=1σxjσxj+1−hN∑j=1σzj, (13)

where are the Pauli matrices at site , is the ferromagnetic exchange coupling, is the transverse external field. Here we assume is even and consider the periodic boundary condition, . This model has a quantum critical point at . Using the Jordan-Wigner transformation,

 σ+j=(σxj+iσyj)2=c†jexp⎛⎝iπ∑l

it can be mapped to a spinless fermion model Dziarmaga2005 (); Franchini2017 (); Dong2018 (),

 H=Nh−2hN∑j=1c†jcj−JN∑j=1(c†j−cj)(c†j+1−cj+1), (15)

where the fermionic operator satisfies either periodic or anti-periodic boundary conditions, , depending on odd or even number of the -quasiparticles, . For simplicity, we confine ourselves to the subspace of even with anti-periodic boundary condition. Using , the fermionic Hamiltonian can be transformed into the bilinear form, , where , , and with .

Below we set to unity and consider the quench protocol where the external field is tuned from its initial value to . The dispersion relation is The characteristic function for transverse Ising chain has an analytical form as shown in Eq. (37) of Appendix B. Figure 2 plots the characteristic function for a quench protocol across the QCP and the corresponding work distribution function obtained through the Gärter-Ellis theorem at different temperatures. As discussed in Appendix B, the work distribution is restricted in the interval, , where is the energy density difference of the highest excited state of the initial phase and the ground state of the final phase and is the energy density difference between the ground state of the initial phase and highest excited state of the final phase. At finite temperature, the small but finite close to and reflects the finite weight of all possible configurations due to thermal effect. However, at zero temperature, since there is no weight for the excited state in the initial phase, is strictly zero if is smaller than the energy density difference of the ground states of the initial and final Hamiltonians, as can be seen in Fig. 2(b). For all the temperatures, we have the normalization condition, , as confirmed in Figs. 2(c) and (d). These results agree with those obtained using the Hamiltonian eigenstate approach Abeling2016 ().

The mean work density and mean irreversible work density at zero temperature can be obtained from Eqs. (7) and (8). Figure 3(a) plots the variation of and as a function of with different lattice size . The mean work density changes continuously with the external field and shows no visible features across the QCP, while a significant peak is seen to grow with increasing in the mean irreversible work density, indicating that the latter is a good indicator of the second order QPT as discussed in Section II. The quantum criticality can then be examined from the scaling behavior in different parameter regimes. In Fig. 3, we see clear logarithmic corrections in depending on , , or , in agreement with the critical exponent of the transverse Ising chain.

We now turn to the double quench protocol where the external field is tuned from to at time and quenched back to at . Figure 4 plots the rate function as a function of for different values of the work density . We see clear cusp in the curve of and continuous variation for other . Such nonanalytic behavior at in the rate function is an indication of the so-called dynamical QPT, which occurs when the time-evolving state is orthogonal to the initial state at certain critical time after quenching a set of control parameters of the Hamiltonian. As discussed in Appendix C for the transverse Ising model, there exists a sequence of critical time, , with

 T∗=π2√1+(h1)2−2h11+h0h1h0+h1. (16)

We should note that such singularity is not present in the the mean (irreversible) work density. As shown in the inset, the mean work density varies smoothly with and is insensitive to the dynamical QPT revealed in . The above results are in good agreement with previous studies using the Loschmidt echo Heyl2013 (), which is defined as and represents the probability amplitude to recover the initial state after the time evolution . For double quench process, the work probability function also gives the return probability to the initial state and could therefore provide a signature when a dynamical QPT occurs Abeling2016 (). The excellent agreement between our results and previous Hamiltonian eigenstate method confirms the validity of our functional field integral approach for both the single and double quench protocols. We may extend it to other models and examine there the effect of quantum criticality on nonequilibrium dynamics and the possible existence of dynamical QPTs under more general circumstances.

### iii.2 The SSH model

In this section, we study the SSH model which was originally proposed for electronic transport in polyacetylene with spontaneous dimerization Su1979 (). Despite its simplicity, the SSH model exhibits a variety of exotic phenomena, such as topological soliton excitation, fractional charge and nontrivial edge states, and has attracted extensive interest in past decades Li2014 (); Shen2012 (); Bernevig2013 (); Asboth2016 (); He2016 (). The model Hamiltonian can be written as

 H=N∑j=1(vc†A,jcB,j+v′c†A,j+1cB,j+h.c.), (17)

where and denote the two sublattices and is the number of unit cells. For even and open boundary condition, there exist two edge modes for but no edge mode for . Thus the model undergoes a topological quantum phase transition at . For periodic boundary condition, one may apply the Fourier transformation, , where or , with . Then the bulk Hamiltonian gets the bilinear form, where and .

Details on the calculations of the mean work density and the mean irreversible work density can be found in the Appendix. Figure 5(a) plots the results as a function of for a single quench from to at . In contrast to that of the transverse Ising model shown in Fig. 3, the mean work density exhibits a clear discontinuity at the critical point. Correspondingly, one observes a sharp resonance-like peak in the mean irreversible work density on a smooth background. As increases, the background evolves into a broad peak, resembling that in the transverse Ising model, but the sharp resonance at the critical point becomes weakened. Except for a small region of the sharp peak around critical , the mean irreversible work density in all other parameter ranges from the relatively smooth background exhibits similar logarithmic scaling with respect to , or . This is plotted in Figs. 5(b-e) and indicates that the critical exponents are , as in the transverse Ising model. In contrast, as one may see from Figs. 6(a) and (b), the mean irreversible work density deviates from the expected scaling at () due to the presence of the sharp resonance. As a matter of fact, we find when is the largest length scale, the scaling relation becomes

 ⟨wirr⟩δ2 ∼lnδ+2N1δ−2λN1δ2,λ−ν>N>δ−ν ∼lnN+2(δ−λ)δ21N,λ−ν>δ−ν>N. (18)

Such non-universal contributions are associated with the topological nature of the SSH model. A straightforward analysis (Appendix C) suggests that they originate from the band crossing at at the topological quantum phase transition. One may separate the contribution from this peculiar point and obtain,

 ⟨wirr⟩k=πδ2=2(δ−λ)δ21N. (19)

For and , in particular, this term has large coefficient and overwhelms the term as shown in Fig. 6(a). To see this more clearly, we plot the mean irreversible work density from momenta and in Figs. 6(c) and (d), respectively. Indeed, there exists a clear and small logarithmic contribution, which is of the same order of magnitude as in Fig. 3(b), but a huge contribution. Intuitively, when the Hamiltonian is quenched across the topological QCP, the necessary existence of the band crossing yields this term. For open boundary condition, these correspond to the edge modes whose contribution to the mean irreversible work density scales inversely with the lattice size .

The result for the double quench protocol is plotted in Fig. 7. Once again, we see a clear cusp in the rate function at at zero temperature, indicating the existence of a dynamical QPT, while the mean work density remains a smooth function of (see the inset). The critical time is now given by with

 T∗=π√1+(v1)2−2v11+v0v1v0+v1, (20)

where and are the values of the Hamiltonian parameter before and after the first quench. We note that since the dynamical QPT in the SSH model is associated with the topological QPT in equilibrium condition, it might be intriguing to think if the nonequilibrium phase might also consist in some topological properties, for example, dynamic edge modes under open boundary condition Y-Wang2017PRB (); Y-Wang2018PRE ().

### iii.3 The BCS model

Now we extend our approach to the superconductivity. Since the model is generally insoluble, we consider here the mean-field BCS Hamiltonian with a bilinear form Altland2010 (),

 H=∑kΨ†k[ξk−Δ−Δ−ξk]Ψk, (21)

where , with the chemical potential , the mean-field order parameter and the Nambu spinor . The parameters and denote the nearest-neighbor and next-nearest-neighbor hoppings on a two-dimensional lattice. We have and with , where is the lattice size along both and directions. Hereafter we set , and .

For single quench, we tune the mean-field parameter from an initial value to such that and . The irreversible work can then be calculated analytically (Appendix C). Figure 8(a) plots as a function of for different lattice size . For large , a logarithmic divergence appears as approaches zero. In this thermodynamical limit (), as shown in Figs. 8(b) and 8(c), we also obtain logarithmic scaling with respect to and . This indicates that the critical exponents are and with .

For the double quench protocol, we change the order parameter from to at and back to at . Using the Gärter-Ellis theorem, we obtain the rate function in Fig. 9 with the initial order parameter . We see a clear nonanalytic behavior in at the critical times , where . Once again, this indicates the existence of a dynamical QPT under double quench. Such a transition only occurs for but is absent for any finite . It must be associated with the superconducting instability. We therefore speculate that the external driven field induces Cooper pair excitations around the initially free electron Fermi surface at . While the calculation of nonequilibrium dynamics might be otherwise involving for a general correlated Hamiltonian, extension of our approach to other strongly correlated phenomena is straightforward under the mean-field approximation.

## Iv Discussion and Conclusions

Our proposal of the functional field integral approach provides an alternate way to calculate the quantum work in an arbitrary time evolution protocol with a general bilinear Hamiltonian. The characteristic function and its cumulants contain all the major information for evaluating the work distribution via the Gärter-Ellis theorem and the mean (irreversible) work density, as well as the fidelity and Loschmidt echo in the double quench process. The applications of our approach to the transverse field Ising model, the SSH model, and the BCS model provide unambiguous evidences for signatures of quantum phase transitions and quantum criticality in the nonequilbirium process. In the sudden quench protocol, this is reflected in the quantum critical scaling of the mean irreversible work density, while in the double quench protocol, dynamical quantum phase transitions may be encoded as a singularity in the rate function. In the SSH model, we also see anomalous corrections due to the topological nature of the quantum phase transition. Compared to the Hamiltonian approach, the functional field integral formalism has the advantage to avoid tedious calculations of the second-order time differential equation and replace them by matrix products at different time intervals, which may be computed efficiently using optimized algorithms at the cost of explicit time evolution of the many-body wave function. Our approach may be extended to general interacting many-body systems in combination with the auxiliary field method to be solved under saddle-point or mean-field approximations. This was shown in our treatment of the BCS Hamiltonian, where a dynamical quantum phase transition is observed in a double quench process. It may be interesting to explore in the future how this dynamical transition is associated with the formation of Copper pairs near the Fermi surface. We also note that only global quench was considered throughout this work, so that translational invariance was kept and the momentum remains a good quantum number. In principle, our approach may also be applied to local quench or random systems by treating the matrix products in real space. An interesting future topic along this line would be the study of many-body localization within the framework of work statistics.

## Acknowledgements

We thank H. T. Quan for the discussions. This work was supported by the National Natural Science Foundation of China (NSFC Grant Nos. 11774401, 11522435), the National Key R&D Program of China (Grant No. 2017YFA0303103), the State Key Development Program for Basic Research of China (Grant No. 2015CB921303), the Strategic Priority Research Program (B) of the Chinese Academy of Sciences (Grant No. XDB07020200) and the Youth Innovation Promotion Association of CAS.

## Appendix A Functional field integral approach to the characteristic function

The Hamiltonian we consider here is

 H(t)=∑kΨ†kAk(t)Ψk, (22)

where and is the vector of the Pauli matrices. Due to the Dirac delta in the work distribution function in Eq. (1), instead of calculating , it is often easier first to calculate its Fourier transformation, namely the characteristic function, . Because of translational invariance, we have with

 Gk(u)=Tr[Uk(0,T0)eiuHk(T0)Uk(T0,0)e−(iu+β)Hk(0)]Tr[e−βHk(0)], (23)

and

 Uk(T0,0)=limM→∞e−iHk(tM,M+1)Δt⋯e−iHk(t1,2)Δt, (24)

where is divided into slices and is an infinitesimal time step as illustrated in Fig. 1. One may then write down the characteristic function by inserting a series of overcomplete bases of the fermionic coherent states on a closed time contour ()Kamenev2011 (),

 ˆ1=∫d[¯¯¯¯ψj,ψj]e−¯¯¯ψjψj∣∣ψj⟩⟨ψj∣∣, (25)

where is defined at with on the contour. This yields

 Gk(u) =1Tr[e−βHk(0)]∫(2M+2∏n=1d[¯¯¯¯ψnk,ψnk])exp(−2M+2∑j=1¯¯¯¯ψjkψjk) ×⟨−ψ2M+2k∣∣eiHk(t2M+1,2M+2)Δt∣∣ψ2M+1k⟩⋯⟨ψM+3k∣∣eiHk(tM+2,M+3)Δt∣∣ψM+2k⟩⟨ψM+2k∣∣eiuHk(T0)∣∣ψM+1k⟩ (26)

Using the definition , we have for any ,

 ⟨ψmk∣∣eαHk(t)Δt∣∣ψnk⟩≈exp(¯¯¯¯ψmkeαΔtAk(t)ψnk). (27)

Thus the partition function at a single model may be reformulated as

 Zk(0) =∫d[¯¯¯¯ψ,ψ]exp[−¯¯¯¯ψ(I+e−βAk(0))ψ]=det(I+e−βAk(0)) =2+Tr[e−βAk(0)], (28)

and similarly,

 Gk(u)=det[I+eiΔtAk(t1,2)⋯eiΔtAk(tM,M+1)eiuAk(T0)e−iΔtAk(tM,M+1)⋯e−iΔtAk(t1,2)e−(iu+β)Ak(0)]det[I+e−βAk(0)], (29)

where is the 22 identity matrix. In deriving above equations, we have used the relation for any 22 matrix and the integral

 ∫d[¯¯¯¯ψ,ψ]exp(−¯¯¯¯ψAψ+¯¯¯¯χψ+¯¯¯¯ψχ)=det(A)exp(¯¯¯¯χA−1χ). (30)

The above formula can be further simplified by defining , with for arbitrary time dependence of between 0 and . We have eventually

 G(u)=∏kdet[I+Bk(T0)]det[I+e−βAk(0)]=∏k2+Tr[Bk(T0)]2+Tr[e−βAk(0)]. (31)

The mean work density is given by the first cumulant of the characteristic function, , yielding

 ⟨w⟩=1Nd∑kTr[(Dk(T0)−Ak(0))e−βAk(0)]2+Tr[e−βAk(0)], (32)

where . On the other hand, taking , one immediately derives the Jarzynski equality Jarzynski1997PRL (); Jarzynski1997PRE (),

 (33)

where is the free energy difference between the final and initial equilibrium states at the inverse temperature , and are the partition functions at and , respectively. Recently, the quantum Jarzynski equality has been experimental verified in trapped ion system An2015 (); Xiong2018 ().

## Appendix B The work statistics and dynamical quantum phase transition

For a single quench where the model Hamiltonian changes from at to at , the characteristic function is

 G(u)=∏k2+Tr[eiuA1ke−(iu+β)A0k]2+Tr[e−βA0k]. (34)

An arbitrary matrix with may be diagonalized under the unitary transformation, , where

 Uk=[μk−ν∗kνkμk],Dk=[Ek00−Ek], (35)

with , , and . The matrix exponential is then

 eαAk=UkeαDkU−1k=⎡⎢ ⎢⎣Ekcosh(αEk)+zksinh(αEk)Ek(xk−iyk)sinh(αEk)Ek(xk+iyk)sinh(αEk)EkEkcosh(αEk)−zksinh(αEk)Ek⎤⎥ ⎥⎦. (36)

We have

 (37)

where .

From Eq. (28), the free energy density difference between the final and initial states in a quench process can be written as,

 (38)

From Eq. (32) and considering for the quench protocol, we have the mean work density,

 ⟨w⟩=1Nd∑kTr[(A1k−A0k)e−βA0k]2+Tr[e−βA0k]=1Nd∑k[(E0k)2−d0k⋅d1k]tanh(12βE0k)E0k. (39)

Combing above equations gives the mean irreversible work density, . For a small quench, , we have the expansion,

 E1k=E0k+x0kδE0k+[(y0k)2+(z0k)2]δ22(E0k)3/2+O(δ3). (40)

Thus at zero temperature

 ⟨w⟩=1Nd∑k−x0kδE0k=−1Nd∑kδ∂Ek∂xk∣∣∣xk=x0k, (41)
 ⟨wirr⟩=1Nd∑k(−x0kδE0k+E1k−E0k)=1Nd∑kδ22∂2Ek∂x2k∣∣ ∣∣xk=x0k. (42)

We see that the singularity in and reflects the first and second-order phase transitions, respectively.

The work distribution may be calculated using the Fourier transformation of . In general, the energy density difference of the ground states, gives the minimal work that can be measured at zero temperature, i.e., for . While at finite temperatures, even the highest excited state may have a small but finite weight due to thermal effect. Thus the minimal work that can be achieved is given by the energy density difference between the ground state of the postquench Hamiltonian with energy density and the highest excited state of the initial Hamiltonian with the energy density , yielding . Oppositely, the maximal work