Functional field integral approach to quantum work

TL;DR

This paper develops a functional field integral method to analyze quantum work statistics in nonequilibrium quantum systems, successfully applying it to models like the transverse Ising chain, SSH model, and BCS superconductivity.

Contribution

It introduces a novel, numerically simple formalism for quantum work analysis applicable to various models, including topological and superconducting systems.

Findings

Correct quantum critical scaling in the transverse Ising chain.

Observation of nonuniversal scaling in the SSH model.

Detection of dynamical quantum phase transitions in multiple protocols.

Abstract

We introduce the functional field integral approach to study the statistics of quantum work under nonequilibrium conditions and derive the general formalism for a bilinear Hamiltonian with arbitrary time dependence. The method is then examined in three models. For the transverse Ising chain, it yields the correct quantum critical scaling and dynamical quantum phase transitions for single and double quench protocols, respectively. For the Su-Schrieffer-Heeger (SSH) model, we observe nonuniversal quantum critical scaling with anomalous $1/N$-correction due to its topological nature. Dynamical quantum phase transitions are observed for three different time evolution protocols but their time periodicity only appears in the double quench case. We then extend our method to the Bardeen-Cooper-Schrieffer (BCS) model for superconductivity and discuss the possibility of its application for general correlated models in combination with either the mean-field approximation or exact Monte Carlo simulations on classical (auxiliary) fields or disorders. Our method has the advantage of numerical simplicity, in the cost of explicit state evolution, and provides a promising way for exploring the physics of quantum work under general conditions.

Equations

$$p(w)=\sum _{n,m}\delta (W-E_m^f+E_n^i)P\left(m^f\mid n^i\right)P\left(n^i\right),$$

$$H\left(t\right)=\sum _k{\Psi }_k^{†}{A}_k\left(t\right){\Psi }_k,$$

$$G(u)={\int }_{-\infty }^{\infty }dw\exp \left(iuw{N}^d\right)p(w)$$

$$=Z_0^{-1}\mathrm{Tr}\left[U^{†}\left(T_0,0\right)e^{iuH\left(T_0\right)}U\left(T_0,0\right)e^{-\left(iu+\beta \right)H\left(0\right)}\right],$$

$$G(u)=\prod _k\frac{2+\mathrm{Tr}\left[B_k\left(T_0\right)\right]}{2+\mathrm{Tr}\left[e^{-\beta A_k\left(0\right)}\right]},$$

$$B_k(T_0)=C_k^{†}\left(T_0\right)e^{iu{A}_k\left(T_0\right)}{C}_k\left(T_0\right)e^{-\left(iu+\beta \right)A_k\left(0\right)},$$

$$⟨w⟩=\frac{1}{N^d}\sum _k\frac{\mathrm{Tr}\left[\left(D_k\left(T_0\right)-A_k\left(0\right)\right)e^{-\beta A_k\left(0\right)}\right]}{2+\mathrm{Tr}\left[e^{-\beta A_k\left(0\right)}\right]},$$

$$\Delta f=-\frac{1}{\beta N^d}\sum _k\ln \frac{{\cosh }^2\left(\beta E_k^1/2\right)}{{\cosh }^2\left(\beta E_k^0/2\right)},$$

$$⟨w⟩=\frac{1}{N^d}\sum _k\frac{\left[{\left(E_k^0\right)}^2-{\mathbf{d}}_k^0\cdot {\mathbf{d}}_k^1\right]\tanh \left(\frac{1}{2}\beta E_k^0\right)}{E_k^0},$$

$$⟨w_{\text{irr}}⟩/{\delta }^2$$

$$\sim N^{2/\nu -\left(d+z\right)},{\delta }^{-\nu }>{\lambda }^{-\nu }>N,$$

$$⟨w_{\text{irr}}⟩/{\delta }^2$$

$$\sim {\lambda }^{\nu \left(d+z\right)-2},N>{\delta }^{-\nu }>{\lambda }^{-\nu }.$$

$$⟨w_{\text{irr}}⟩/{\delta }^2$$

$$\sim N^{2/\nu -\left(d+z\right)},{\lambda }^{-\nu }>{\delta }^{-\nu }>N.$$

$$r(w)=-\underset{R\in \mathbb{R}}{\inf }[wR-c(R)],$$

$$C_k\left(T_0\right)\approx e^{-i\Delta t{A}_k\left(\frac{t_{M+1}+t_M}{2}\right)}\cdots e^{-i\Delta t{A}_k\left(\frac{t_1+t_2}{2}\right)},$$

$$H=-J\sum _{j=1}^N{\sigma }_j^x{\sigma }_{j+1}^x-h\sum _{j=1}^N{\sigma }_j^z,$$

$${\sigma }_j^{+}=\frac{\left({\sigma }_j^x+i{\sigma }_j^y\right)}{2}=c_j^{†}\exp \left(i\pi \sum _{l<j}{c}_l^{†}{c}_l\right),$$

$$H=Nh-2h\sum _{j=1}^N{c}_j^{†}{c}_j-J\sum _{j=1}^N\left(c_j^{†}-c_j\right)\left(c_{j+1}^{†}-c_{j+1}\right),$$

$$T^{\ast }=\frac{\pi }{2\sqrt{1+{\left(h_1\right)}^2-2h_1\frac{1+h_0{h}_1}{h_0+h_1}}}.$$

$$H=\sum _{j=1}^N\left(v{c}_{A,j}^{†}{c}_{B,j}+v^{\prime }{c}_{A,j+1}^{†}{c}_{B,j}+h.c.\right),$$

$$\frac{⟨w_{\text{irr}}⟩}{{\delta }^2}$$

$$\sim \ln N+\frac{2\left(\delta -\lambda \right)}{{\delta }^2}\frac{1}{N},{\lambda }^{-\nu }>{\delta }^{-\nu }>N.$$

$$\frac{{⟨w_{\text{irr}}⟩}_{k=\pi }}{{\delta }^2}=\frac{2\left(\delta -\lambda \right)}{{\delta }^2}\frac{1}{N}.$$

$$T^{\ast }=\frac{\pi }{\sqrt{1+{\left(v_1\right)}^2-2v_1\frac{1+v_0{v}_1}{v_0+v_1}}}.$$

$$H=\sum_{k}\Psi_{k}^{\dagger}\left[\begin{array}[c]{cc}\xi_{k}&-\Delta\\ -\Delta&-\xi_{k}\end{array}\right]\Psi_{k},$$

$$G(u)=\frac{{\sum }_{\{S_j^z\}}\det \left[I+B(\{S_j^z\},T_0)\right]}{{\sum }_{\{S_j^z\}}\det \left[I+B(\{S_j^z\},0)\right]},$$

$$H\left(t\right)=\sum _k{\Psi }_k^{†}{A}_k\left(t\right){\Psi }_k,$$

$$G_k(u)=\frac{\mathrm{Tr}\left[U_k\left(0,T_0\right)e^{iu{H}_k\left(T_0\right)}{U}_k\left(T_0,0\right)e^{-\left(iu+\beta \right)H_k\left(0\right)}\right]}{\mathrm{Tr}\left[e^{-\beta H_k\left(0\right)}\right]},$$

$$U_k\left(T_0,0\right)=\underset{M\to \infty }{\lim }{e}^{-i{H}_k\left(t_{M,M+1}\right)\Delta t}\cdots e^{-i{H}_k\left(t_{1,2}\right)\Delta t},$$

$$\hat{1}=\int d\left[{\overline{\psi }}^j,{\psi }^j\right]e^{-{\overline{\psi }}^j{\psi }^j}|{\psi }^j⟩⟨{\psi }^j|,$$

$$G_k\left(u\right)$$

Full text

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

[TABLE]

where $w=W/N^{d}$ is the work density, $N$ is the lattice size, $d$ is the dimensionality, $E_{n}^{i}$ is the $n$-th eigenvalue of the initial state ($|n^{i}\rangle$), and $E_{m}^{f}$ is the $m$-th eigenvalue of the final state ($|m^{f}\rangle$). $P(n^{i})=\exp(-\beta E_{n}^{i})/Z(0)$ denotes the probability distribution of the initial canonical state, and $P(m^{f}\mid n^{i})\equiv|\langle m^{f}|U(T_{0},0)|n^{i}\rangle|^{2}$ accounts for the transition probability between the initial and final states governed by the time-dependent evolution, $U(T_{0},0)=\mathcal{T}\exp\left[-\operatorname*{i}\int_{0}^{T_{0}}dtH(t)\right]$, where $\mathcal{T}$ is the time ordering operator and $T_{0}$ is the period for the time evolution. The reduced Planck constant $\hbar$ is set to unity. The nonequilibrium dynamics is fully incorporated in the time-dependent Hamiltonian $H(t)$ 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 hard to be extended to arbitrary time dependence and 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, which is a many-body extension of the path integral approach to quantum fields Altland2010 . Our approach 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,

[TABLE]

where $A_{k}(t)=\mathbf{d}_{k}(t)\cdot\boldsymbol{\sigma}$ is a matrix and $\boldsymbol{\sigma}$ 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 :

[TABLE]

where $Z_{0}=\operatorname*{Tr}\operatorname{e}^{-\beta H\left(0\right)}$ is the partition function at $t=0$. The time contour of the path integral in $G(u)$ is illustrated in Fig. 1. After some tedious calculations using the functional field integral techniques (see Appendix A), we obtain

[TABLE]

where

[TABLE]

with $C_{k}\left(T_{0}\right)=\mathcal{T}\exp\left[-\operatorname*{i}\int_{0}^{T_{0}}dtA_{k}\left(t\right)\right]$. The mean work density is given by the first cumulant of the characteristic function, $\left\langle w\right\rangle=-\operatorname*{i}dG(u)/(N^{d}du)|_{u=0}$, yielding

[TABLE]

with $D_{k}(T_{0})=C_{k}^{\dagger}(T_{0})A_{k}(T_{0})C_{k}(T_{0})$, where $C_{k}(T_{0})$ can be computed numerically for arbitrary time dependence.

Note that the mean work density $\left\langle w\right\rangle$ always exceeds the free energy density difference $\Delta f$ between the initial and final equilibrium states (both with the same inverse temperature $\beta$), as stated in the second law of thermodynamics. A mean irreversible work density can thus be defined as their difference, $\langle w_{\text{irr}}\rangle=\langle w\rangle-\Delta f\geq 0$, which characterizes the irreversibility of the nonequilibrium process and is directly related to the entropy increase between the final and initial equilibrium states, $\Delta s=\beta\langle w_{\text{irr}}\rangle$, for a closed quantum system without heat transfer Shraddha2015 .

The mean work density $\left\langle w\right\rangle$ and the mean irreversible work density $\left\langle w_{\text{irr}}\right\rangle$ can be used to identify equilibrium phase transitions. Considering a quench process where the model Hamiltonian changes from $A_{k}^{0}=\mathbf{d}_{k}^{0}\cdot\boldsymbol{\sigma}$ at time $t=0^{-}$ to $A_{k}^{1}=\mathbf{d}_{k}^{1}\cdot\boldsymbol{\sigma}$ at $t=T_{0}=0^{+}$, we have $D_{k}(0^{+})=A_{k}^{1}$. As shown in Appendix B, one can derive an explicit formula for the free energy density difference,

[TABLE]

and the mean work density,

[TABLE]

where $E_{k}^{0,1}=\left|\mathbf{d}_{k}^{0,1}\right|$. 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, $\mathbf{d}_{k}^{0}=\left(x_{k}^{0},y_{k}^{0},z_{k}^{0}\right)\rightarrow\mathbf{d}_{k}^{1}=\left(x_{k}^{0}+\delta,y_{k}^{0},z_{k}^{0}\right)$. It can be shown analytically (Appendix B) that $\left\langle w\right\rangle=-\sum_{k}\frac{\delta}{N^{d}}\frac{\partial E_{k}}{\partial x_{k}}|_{x_{k}=x_{k}^{0}}$ and $\left\langle w_{\text{irr}}\right\rangle=\sum_{k}\frac{\delta^{2}}{2N^{d}}\frac{\partial^{2}E_{k}}{\partial x_{k}^{2}}|_{x_{k}=x_{k}^{0}}$ 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 $\left\langle w\right\rangle$ and $\left\langle w_{\text{irr}}\right\rangle$ reflect the first or second-order phase transitions.

For a (second-order) QPT, scaling analysis of $\left\langle w_{\text{irr}}\right\rangle$ can yield key information on the critical exponents Shraddha2015 . In the heat susceptibility limit, where $\delta^{-v}$ is the largest length scale, the mean irreversible work density at zero temperature scales as

[TABLE]

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

[TABLE]

On the other hand, if the system is prepared very close to the QCP such that $\lambda^{-\nu}$ is larger than all other length scales, the scaling relation becomes

[TABLE]

Note that when the combination $\nu\left(d+z\right)$ exceeds 2, the scaling of $\left\langle w_{\text{irr}}\right\rangle$ is non-universal. In the marginal case, $\nu\left(d+z\right)=2$, 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 $A_{k}^{0}$ to $A_{k}^{1}$ at $t=0$ and then back to $A_{k}^{0}$ for $t\geq T_{0}$, can be used to describe dynamical QPTs Heyl2013 . In this case, the free energy density difference $\Delta f$ is zero and we have $\langle w_{\text{irr}}\rangle=\langle w\rangle$ and $D_{k}(T_{0})=e^{\operatorname*{i}T_{0}A_{k}^{1}}A_{k}^{0}e^{-\operatorname*{i}T_{0}A_{k}^{1}}$. 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 $G(u)$, 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, $p(w)\sim e^{-N^{d}r(w)}$, which defines the rate function $r(w)\geq 0$. The theorem states that $r(w)$ can be obtained via a Legendre-Fenchel transformation,

[TABLE]

where $c(R)=-\lim_{N\rightarrow\infty}N^{-d}\ln G(u=\operatorname*{i}R)$ is a scaled cumulant generating function assumed to be differentiable with respect to the real variable $R\in\mathbb{R}$, and the infimum is evaluated within the domain of definition of $c(R)$ including $R=\pm\infty$. The dynamical QPT is then manifested as a singular point of $r(w)$, as discussion in Appendix B.

For general time dependence, $C_{k}(T_{0})$ may be computed numerically using the matrix products,

[TABLE]

which becomes accurate as $M\rightarrow\infty$ and $\Delta t=T_{0}/M\rightarrow 0$. If all $A_{k}(t)$ commute, the above formula reduces to a simple time integral of $A_{k}(t)$ in the exponent. It is anticipated that detailed investigations for arbitrary time dependence might reveal more interesting properties of nonequilibrium physics. In this work, we will focus on quench protocols for the sake of simplicity and only give as an example a brief discussion on general cases in Section III.2.

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

[TABLE]

where $\sigma_{j}^{\alpha}\,(\alpha=x,z)$ are the Pauli matrices at site $j$, $J>0$ is the ferromagnetic exchange coupling, $h$ is the transverse external field. Here we assume $N$ is even and consider the periodic boundary condition, $\sigma_{N+1}^{\alpha}=\sigma_{1}^{\alpha}$. This model has a quantum critical point at $h=J$. Using the Jordan-Wigner transformation,

[TABLE]

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

[TABLE]

where the fermionic operator $c_{j}$ satisfies either periodic or anti-periodic boundary conditions, $c_{N+1}=\pm c_{1}$, depending on odd or even number of the $c$-quasiparticles, $M=\sum_{j=1}^{N}c_{j}^{\dagger}c_{j}$. For simplicity, we confine ourselves to the subspace of even $M$ with anti-periodic boundary condition. Using $c_{q}=N^{-1/2}\sum_{j=1}^{N}c_{j}\exp(\operatorname*{i}qj)$, the fermionic Hamiltonian can be transformed into the bilinear form, $H=\sum_{k>0}\Psi_{k}^{\dagger}\mathbf{d}_{k}\cdot\boldsymbol{\sigma}\Psi_{k}$, where $\Psi_{k}^{\dagger}=(c_{k}^{\dagger},c_{-k})$, $\mathbf{d}_{k}=(0,-2J\sin k,-2h-2J\cos k)$, and $k=\pm\pi(2m-1)/N$ with $m=1,\ldots,\,N/2$.

Below we set $J$ to unity and consider the quench protocol where the external field $h$ is tuned from its initial value $h_{0}$ to $h_{1}=h_{0}+\delta$. The dispersion relation is $\epsilon_{k}(h)=|\mathbf{d}_{k}|=2\sqrt{1+h^{2}+2h\cos k}.$ The characteristic function $G(u)=\prod_{k>0}G_{k}(u)$ for transverse Ising chain has an analytical form as shown in Eq. (39) 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, $[w_{\min},w_{\max}]$, where $w_{\min}$ is the energy density difference of the highest excited state of the initial phase and the ground state of the final phase and $w_{\max}$ 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 $p(w)$ close to $w_{\min}$ and $w_{\max}$ 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, $p(w)$ is strictly zero if $w$ 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, $G(0)=\int_{-\infty}^{\infty}dw\,p(w)=1$, 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 $\langle w\rangle/\delta$ and $\langle w_{\text{irr}}\rangle/\delta^{2}$ as a function of $h_{0}$ with different lattice size $N$. The mean work density changes continuously with the external field $h_{0}$ and shows no visible features across the QCP, while a significant peak is seen to grow with increasing $N$ 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 $\langle w_{\text{irr}}\rangle$ depending on $N$, $\lambda$, or $\delta$, in agreement with the critical exponent $\nu=z=1$ of the transverse Ising chain.

We now turn to the double quench protocol where the external field $h$ is tuned from $h_{0}<1$ to $h_{1}>1$ at time $t=0$ and quenched back to $h_{0}$ at $t=T_{0}$. Figure 4 plots the rate function as a function of $T_{0}$ for different values of the work density $w=W/N$. We see clear cusp in the curve of $w=0$ and continuous variation for other $w$. Such nonanalytic behavior at $w=0$ in the rate function is an indication of the so-called dynamical QPT, which occurs when the time-evolving state $|\psi(t)\rangle$ 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, $t_{c}=(n+1/2)T^{\ast}$, with

[TABLE]

Its inverse, $\omega^{\ast}=2\pi/T^{\ast}$, seems to define a characteristic excitation energy for the transition Piccitto2019 . We should note that such singularity is not present in the mean (irreversible) work density. As shown in the inset, the mean work density varies smoothly with $T_{0}$ and is insensitive to the dynamical QPT revealed in $r(w,T_{0})$. The above results are in good agreement with previous studies using the Loschmidt echo Heyl2013 , which is defined as $L\left(t\right)=\left|\left\langle\psi\left|U\left(t\right)\right|\psi\right\rangle\right|^{2}$ and represents the probability amplitude to recover the initial state after the time evolution $U\left(t\right)$. For double quench process, the work probability function $p\left(w=0,T_{0}\right)$ 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 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

[TABLE]

where $A$ and $B$ denote the two sublattices and $N$ is the number of unit cells. For even $N$ and open boundary condition, there exist two edge modes for $v<v^{\prime}$ but no edge mode for $v>v^{\prime}$. Thus the model undergoes a topological quantum phase transition at $v=v^{\prime}$. For periodic boundary condition, one may apply the Fourier transformation, $c_{\alpha,j}=N^{-1/2}\sum_{k}c_{\alpha,k}\exp\left(\operatorname*{i}kj\right)$, where $\alpha=A$ or $B$, $k=2m\pi/N$ with $m=1-N/2,\ldots,\,N/2$. Then the bulk Hamiltonian gets the bilinear form, $H=\sum_{k}\Psi_{k}^{\dagger}\mathbf{d}_{k}\cdot\boldsymbol{\sigma}\Psi_{k},$ where $\Psi_{k}^{\dagger}=\left(c_{A,k}^{\dagger},c_{B,k}^{\dagger}\right)$ and $\mathbf{d}_{k}=\left(v+\cos k,\sin k,0\right)$.

Details on the calculations of the mean work density and the mean irreversible work density can be found in the Appendix C. Figure 5(a) plots the results as a function of $v_{0}$ for a single quench from $v=v_{0}$ to $v_{0}+\delta$ at $v^{\prime}=1$. 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 $N$ 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 $v_{0}=1$, the mean irreversible work density in all other parameter ranges from the relatively smooth background exhibits similar logarithmic scaling with respect to $N$, $\lambda$ or $\delta$. This is plotted in Figs. 5(b-e) and indicates that the critical exponents are $\nu=z=1$, 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 $v_{0}=1$ ($\lambda=0$) due to the presence of the sharp resonance. As a matter of fact, we find when $\lambda^{-\nu}$ is the largest length scale, the scaling relation becomes

[TABLE]

Such non-universal $1/N$ 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 $k=\pi$ at the topological quantum phase transition. One may separate the contribution from this peculiar point and obtain,

[TABLE]

For $\lambda^{-\nu}>\delta^{-\nu}>N$ and $\nu=1$, in particular, this $1/N$ term has large coefficient and overwhelms the $\ln N$ term as shown in Fig. 6(a). To see this more clearly, we plot the mean irreversible work density from momenta $k\neq\pi$ and $k=\pi$ 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 $1/N$ contribution. Intuitively, when the Hamiltonian is quenched across the topological QCP, the necessary existence of the band crossing yields this $1/N$ 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 $N$.

The rate function for the double quench protocol of the SSH model is plotted in Fig. 7 and actually similar to that of the Ising model (up to a factor of 2) since both models have the same bilinear Hamiltonian except for the boundary conditions. Thus their difference is delicate and actually diminishes in the thermodynamic limit. For comparison, we also plot the results for two other protocols with different time dependent $v(t)$ as shown in the inset. Interestingly, we see clear cusps in $r(w,T_{0})$ at $w=0$ and zero temperature in all three cases, indicating the possibility of general existence of dynamical QPTs. However, no time periodicity is seen except for the double quench case, where the critical time is given by $t_{c}=\left(n+1/2\right)T^{\ast}$ with

[TABLE]

More elaborate investigations might be able to reveal the true controlling parameter for the dynamical QPT. We note that since the dynamical QPT in the SSH model is associated with the equilibrium QPT with topological properties, 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 conditions 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 and discuss later in section IV the possibility of going beyond the mean-field approximation. The Hamiltonian is Altland2010 ,

[TABLE]

where $\xi_{k}=-2v[\cos(k_{x})+\cos(k_{y})]-4v^{\prime}\cos(k_{x})\cos(k_{y})-\mu$, with the chemical potential $\mu$, the mean-field order parameter $\Delta$ and the Nambu spinor $\Psi_{k}^{\dagger}=(c_{k\uparrow}^{\dagger},c_{-k\downarrow})$. The parameters $v$ and $v^{\prime}$ denote the nearest-neighbor and next-nearest-neighbor hoppings on a two-dimensional lattice. We have $k_{x}=2\pi m/N$ and $k_{y}=2\pi l/N$ with $m,l=1-N/2,\ldots,\,N/2$, where $N$ is the lattice size along both $x$ and $y$ directions. Hereafter we set $v=0.435$, $v^{\prime}=0.05$ and $\mu=0.5$.

For single quench, we tune the mean-field parameter $\Delta$ from an initial value $\Delta_{0}$ to $\Delta_{1}=\Delta_{0}+\delta$ such that $\mathbf{d}_{k}^{0}=(-\Delta_{0},0,\xi_{k})$ and $\mathbf{d}_{k}^{1}=(-\Delta_{1},0,\xi_{k})$. The irreversible work can then be calculated analytically (Appendix C). Figure 8(a) plots $\langle w_{\text{irr}}\rangle/\delta^{2}$ as a function of $\Delta_{0}$ for different lattice size $N$. For large $N$, a logarithmic divergence appears as $\Delta_{0}$ approaches zero. In this thermodynamical limit ($N=\infty$), as shown in Figs. 8(b) and 8(c), we also obtain logarithmic scaling with respect to $\lambda$ and $\delta$. This indicates that the critical exponents are $\nu=1/2$ and $z=2$ with $d=2$.

For the double quench protocol, we change the order parameter from $\Delta_{0}$ to $\Delta_{1}$ at $t=0$ and back to $\Delta_{0}$ at $t=T_{0}$. Using the Gärter-Ellis theorem, we obtain the rate function $r(w,T_{0})$ in Fig. 9 with the initial order parameter $\Delta_{0}=0$. We see a clear nonanalytic behavior in $r(w=0,T_{0})$ at the critical times $t_{c}=(n+1/2)T^{\ast}$, where $T^{\ast}=\pi/\Delta_{1}$. Once again, this indicates the existence of a dynamical QPT under double quench. Such a transition only occurs for $\Delta_{0}=0$ but is absent for any finite $\Delta_{0}$. 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 $\Delta_{0}=0$. 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 $1/N$ 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. As an example, we calculate the rate function for the SSH model with different time evolution protocols and find dynamical quantum phase transitions in all cases. However, the time periodicity only appears in the double quench case. Our approach provides a possibility for revealing the true driving parameter of the dynamical transition under general nonequilibrium conditions.

The transverse Ising model represents the few examples that can be mapped to a simple bilinear form. Often, an interacting many-body system cannot be solved even under equilibrium conditions, and approximations or numerical simulations have to be used with specially designed algorithms. These include the Hubbard model for cuprates and pnictides and the periodic Anderson model or the Kondo lattice model for heavy fermions. In these cases, certain auxiliary fields, based also on the functional field integral formalism, might be introduced to decouple the interacting terms in the Hamiltonian, so that the model can be solved under the mean-field (or saddle-point) approximations or simulated exactly using efficient numerical techniques Gubernatis2016 . The mean-field approximation can often capture some of the essential physics of the correlated model. As shown for the BCS Hamiltonian, a dynamical quantum phase transition is observed in a double quench process, which might be an intrinsic property of the model ascribed to the formation of Copper pairs near the Fermi surface driven by the time dependent pairing force.

To go beyond the mean-field approximation, one may consider Monte Carlo simulations of the auxiliary fields. A simpler situation is that the electrons are coupled to classical fields or subject to disorders but the Hamiltonian remains a bilinear form for each field or disorder configuration, as is in the Anderson localization model or the Ising Kondo lattice model. For the latter, we may have $H(t)=-v\sum_{\langle ij\rangle,\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}+\text{h.c.})+J(t)\sum_{j\sigma}S_{j}^{z}\sigma c_{j\sigma}^{\dagger}c_{j\sigma}$, where $S_{j}^{z}$ is the Ising spin defined at site $j$ and $J(t)$ is the time-dependent coupling. Since all $S_{j}^{z}$ commute with $H$, the spins do not evolve with time. For each spin configuration, the electron Hamiltonian has a bilinear form in real space, $H(\{S_{j}^{z}\},t)=\sum_{i\sigma,l\sigma^{\prime}}c_{i\sigma}^{\dagger}A_{i\sigma,l\sigma^{\prime}}(\{S_{j}^{z}\},t)c_{l\sigma^{\prime}}$. Thus our derivation in appendix A can be extended to give a similar formula after the electron degrees of freedom are integrated out,

[TABLE]

where $B$ is defined as in Eq. (5) but with $A_{k}(t)$ replaced by the matrix $A(\{S_{j}^{z}\},t)$ for each spin configuration $\{S_{j}^{z}\}$. Then the Monte Carlo simulations can be applied on the Ising spin configurations Zhong2019 , with matrix products in real space instead of momentum space due to the lack of translation symmetry. Similar procedure may be extended to more complicated cases such as the Hubbard model or periodic Anderson model, where the local Coulomb interaction may be decoupled by introducing an auxiliary field $s_{i}$ at each site and time slice such that $e^{-\Delta\tau U(n_{i\uparrow}-1/2)(n_{i\downarrow}-1/2)}=c\sum_{s_{i}=\pm 1}e^{\alpha s_{i}(n_{i\uparrow}-n_{i\downarrow})}$, where $c$ and $\alpha$ are constants depending on $\Delta\tau U$. Then one may perform Monte Carlo sampling on the auxiliary field configurations Gubernatis2016 . In both cases, detailed numerical realization may be limited by severe sign problem or other issues and thus require careful design of the algorithms. Nevertheless, our approach seems to provide a general scheme for the exploration of quantum work in nonequilibrium correlated systems. One interesting 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 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), and the National Youth Top-Notch Talent Support Program of China.

Appendix A Functional field integral approach to the characteristic function

The Hamiltonian we consider here is

[TABLE]

where $A_{k}(t)=\mathbf{d}_{k}(t)\cdot\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}$ is the vector of the Pauli matrices. Due to the Dirac delta in the work distribution function in Eq. (1), instead of calculating $p(w)$, it is often easier first to calculate its Fourier transformation, namely the characteristic function, $G(u)=\int_{-\infty}^{\infty}dw\,\text{e}^{\operatorname*{i}uwN^{d}}p(w)$. Because of translational invariance, we have $G(u)=\prod\limits_{k}G_{k}(u)$ with

[TABLE]

and

[TABLE]

where $T_{0}$ is divided into $M$ slices and $\Delta t=T_{0}/M$ 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 ($0\rightarrow T_{0}\rightarrow 0$)Kamenev2011 ,

[TABLE]

where $|\psi^{j}\rangle$ is defined at $t_{j}$ with $j=1,2,\ldots,2M+2$ on the contour. This yields

[TABLE]

Using the definition $\langle\psi_{k}^{m}|\psi_{k}^{n}\rangle=\exp\left(\overline{\psi}_{k}^{m}\psi_{k}^{n}\right)$, we have for any $t$,

[TABLE]

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

[TABLE]

and similarly,

[TABLE]

where $I$ is the 2$\times$2 identity matrix. In deriving above equations, we have used the integral

[TABLE]

The above formula can be further simplified by defining $B_{k}(T_{0})=C_{k}^{\dagger}(T_{0})\operatorname{e}^{\operatorname*{i}uA_{k}(T_{0})}C_{k}(T_{0})\operatorname{e}^{-(\operatorname*{i}u+\beta)A_{k}(0)}$, with $C_{k}(T_{0})=\mathcal{T}\exp\left[-\operatorname*{i}\int_{0}^{T_{0}}dtA_{k}(t)\right]$ for arbitrary time dependence of $A_{k}(t)$ between 0 and $T_{0}$. Note that $B_{k}(0)=\operatorname{e}^{-\beta A_{k}(0)}$. We have eventually

[TABLE]

following the identity, $\det(I+A)=1+\operatorname*{Tr}A+\det A$, for any 2$\times$2 matrix $A$. The mean work density is given by the first cumulant of the characteristic function, $\langle w\rangle=-\operatorname*{i}dG(u)/(N^{d}du)|_{u=0}$, yielding

[TABLE]

where $D_{k}(T_{0})=C_{k}^{\dagger}(T_{0})A_{k}(T_{0})C_{k}(T_{0})$. On the other hand, taking $u=\operatorname*{i}\beta$, one immediately derives the Jarzynski equality Jarzynski1997PRL ; Jarzynski1997PRE ,

[TABLE]

where $\Delta F=F_{1}-F_{0}$ is the free energy difference between the final and initial equilibrium states at the inverse temperature $\beta$, $Z(0)$ and $Z(T_{0})$ are the partition functions at $t=0$ and $t=T_{0}$, 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 $A_{k}^{0}=\mathbf{d}_{k}^{0}\cdot\boldsymbol{\sigma}$ at $t=0^{-}$ to $A_{k}^{1}=\mathbf{d}_{k}^{1}\cdot\boldsymbol{\sigma}$ at $t=T_{0}=0^{+}$, the characteristic function is

[TABLE]

An arbitrary matrix $A_{k}=\mathbf{d}_{k}\cdot\boldsymbol{\sigma}$ with $\mathbf{d}_{k}=(x_{k},y_{k},z_{k})$ may be diagonalized under the unitary transformation, $U_{k}^{-1}A_{k}U_{k}=D_{k}$, where

[TABLE]

with $E_{k}=|\mathbf{d}_{k}|$, $\mu_{k}=\sqrt{(E_{k}+z_{k})/2E_{k}}$, and $\nu_{k}=\frac{x_{k}+\operatorname*{i}y_{k}}{\sqrt{x_{k}^{2}+y_{k}^{2}}}\sqrt{\frac{E_{k}-z_{k}}{2E_{k}}}$. The matrix exponential $\operatorname{e}^{\alpha A_{k}}$ is then

[TABLE]

We have

[TABLE]

where $E_{k}^{0,1}=|\mathbf{d}_{k}^{0,1}|$.

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

[TABLE]

From Eq. (34) and considering $D_{k}(T_{0})=A_{k}^{1}$ for the quench protocol, we have the mean work density,

[TABLE]

Combing above equations gives the mean irreversible work density, $\langle w_{\text{irr}}\rangle=\langle w\rangle-\Delta f$. For a small quench, $\mathbf{d}_{k}^{0}=(x_{k}^{0},y_{k}^{0},z_{k}^{0})$ $\rightarrow$ $\mathbf{d}_{k}^{1}=(x_{k}^{0}+\delta,y_{k}^{0},z_{k}^{0})$, we have the expansion,

[TABLE]

Thus at zero temperature

[TABLE]

[TABLE]

We see that the singularity in $\left\langle w\right\rangle$ and $\left\langle w_{\text{irr}}\right\rangle$ reflects the first and second-order phase transitions, respectively.

The work distribution $p(w)$ may be calculated using the Fourier transformation of $G(u)$. In general, the energy density difference of the ground states, $\Delta_{\min}=\Delta f=-\frac{1}{N^{d}}\sum_{k}\left[E_{k}^{1}-E_{k}^{0}\right]$ gives the minimal work that can be measured at zero temperature, i.e., $p(w<\Delta_{\min})=0$ for $T=0$. While at finite temperatures, even the highest excited state may have a small but finite weight due to thermal effect. Thus the minimal work $w_{\min}$ that can be achieved is given by the energy density difference between the ground state of the postquench Hamiltonian with energy density $\epsilon^{g}=-\frac{1}{N^{d}}\sum_{k}E_{k}^{1}$ and the highest excited state of the initial Hamiltonian with the energy density $\epsilon^{h}=\frac{1}{N^{d}}\sum_{k}E_{k}^{0}$, yielding $w_{\min}=-\frac{1}{N^{d}}\sum_{k}(E_{k}^{1}+E_{k}^{0})$. Oppositely, the maximal work $w_{\max}$ is given by the energy density difference between the highest excited state of the postquench Hamiltonian with energy density $\epsilon^{h}=\frac{1}{N^{d}}\sum_{k}E_{k}^{1}$ and the ground excited state of the initial Hamiltonian with the energy density $\epsilon^{g}=-\frac{1}{N^{d}}\sum_{k}E_{k}^{0}$, yielding $w_{\max}=\frac{1}{N^{d}}\sum_{k}(E_{k}^{1}+E_{k}^{0})=-w_{\min}$. Beyond the interval $[w_{\min},w_{\max}]$, the distribution function $p(w)=0$ and the rate function $r(w)=\infty$ at all temperatures.

For a double quench, where the Hamiltonian changes from $A_{k}^{0}$ to $A_{k}^{1}$ at $t=0$ and back to $A_{k}^{0}$ for $t\geq T_{0}$, we have $C_{k}(T_{0})=\operatorname{e}^{-\operatorname*{i}T_{0}A_{k}^{1}}$. The characteristic function at zero temperature can be evaluated to be

[TABLE]

Because $c(R)=-\lim_{N\rightarrow\infty}\frac{1}{N^{d}}\ln G(u=\operatorname*{i}R)$ and $r(w=0)=\sup_{R\in\mathbb{R}}c(R)$, the singularity in the rate function is associated with the roots of $G(R)$. Thus a dynamical quantum phase transition may occur at the critical time, $t_{c}=\frac{\left(2n+1\right)\pi}{2E_{k_{c}}^{1}}$ ($n=0$, 1, 2, …) if there exits a critical momentum $k_{c}$ satisfying $\mathbf{d}_{k_{c}}^{0}\cdot\mathbf{d}_{k_{c}}^{1}=0$. Under this condition, we have at $T_{0}=t_{c}$,

[TABLE]

Then the singularity occurs when $G(R)=\prod_{k}G_{k}(R)=0$ as $R\rightarrow\infty$. While for all other $T_{0}\neq t_{c}$ or $k\neq k_{c}$, there exists a term in the brace proportional to $\operatorname{e}^{RE_{k_{c}}^{0}}$, which will cancel the prefactor $\operatorname{e}^{-RE_{k}^{0}}$ and produce a nonzero $G(R)$ for all $R$.

Appendix C Application to the models

For the transverse Ising chain, we have $\mathbf{d}_{k}=(0,-2\sin k,-2h-2\cos k)$. Using the dispersion relation $\epsilon_{k}(h)=|\mathbf{d}_{k}|=2\sqrt{1+h^{2}+2h\cos k}$, we can calculate $\left\langle w_{\text{irr}}\right\rangle$ for the single quench from $h_{0}<1$ to $h_{0}>1$,

[TABLE]

Using $t_{c}=\frac{(2n+1)\pi}{2E_{k_{c}}^{1}}$ and $\mathbf{d}_{k_{c}}^{0}\cdot\mathbf{d}_{k_{c}}^{1}=0$ for the value of $k_{c}$, the critical times for the dynamical QPT during a double quench process are

[TABLE]

For the SSH model, we have $\mathbf{d}_{k}=\left(v+\cos k,\sin k,0\right)$ and $\epsilon_{k}\left(v\right)=$ $\left|\mathbf{d}_{k}\right|=\sqrt{1+v^{2}+2v\cos k}$. Thus for a single quench from $v_{0}<1$ to $v_{1}=v_{0}+\delta>1$,

[TABLE]

Because of the topological nature of the QPT in the SSH model, there exists a band crossing at $k=\pi$ at the QCP. We can thus separate the mean irreversible work density into two terms, $\langle w_{\text{irr}}\rangle=\langle w_{\text{irr}}\rangle_{k\neq\pi}+\left\langle w_{\text{irr}}\right\rangle_{k=\pi}$. The results are shown in Fig. 6. While the first term gives the usual $\ln N$ scaling due to quantum criticality, the second term $\langle w_{\text{irr}}\rangle_{k=\pi}=\frac{2}{N}\left(\delta-\lambda\right)$ yields an additional $1/N$ contribution, which becomes dominant and leads to the anomalous $1/N$ scaling when $\lambda^{-\nu}$ is the largest length scale. For the double quench process, the critical times for the dynamical QPT are

[TABLE]

For the BCS model with $\mathbf{d}_{k}=(-\Delta,0,\xi_{k})$, the mean irreversible work density for the single quench process is

[TABLE]

and the dynamical QPT during a double quench process with $\Delta_{0}=0$ occurs at

[TABLE]