Work statistics for Quantum Spin Chains: characterizing quantum phase transitions, benchmarking time evolution, and examining passivity of quantum states

We study three aspects of work statistics in the context of the fluctuation theorem for the quantum spin chains up to $1024$ sites by numerical methods based on matrix-product states (MPS). First, we use our numerical method to evaluate the moments/cumulants of work done by sudden quench process on the Ising or Haldane spin chains and study their behaviors across the quantum phase transitions. Our results show that, up to the fourth cumulant, the work statistics can indicate the quantum phase transition characterized by the local order parameters but barely for purely topological phase transitions. Second, we propose to use the fluctuation theorem, such as Jarzynski's equality, which relates the real-time correlator to the ratio of the thermal partition functions, as a benchmark indicator for the numerical real-time evolving methods. Third, we study the passivity of ground and thermal states of quantum spin chains under some cyclic impulse processes. We show that the passivity of thermal states and ground states under the hermitian actions are ensured by the second laws and variational principles, respectively, and also verify it by numerical calculations. Besides, we also consider the passivity of ground states under non-hermitian actions, for which the variational principle cannot be applied. Despite that, we find no violation of passivity from our numerical results for all the cases considered in the Ising and Haldane chains. {Overall, we demonstrate that the work statistics for the sudden quench and impulse processes can be evaluated precisely by the numerical MPS method to characterize quantum phase transitions and examine the passivity of quantum states. We also propose to exploit the universality of the fluctuation theorem to benchmark the numerical real-time evolutions in an algorithm and model independent way.

IV. Numerical results for characterizing the quantum phase transitions 8 A. Quantum spin-1/2 chain in magnetic field 8 B. Anisotropic quantum spin-1 chain 9 V. Benchmarking the real-time evolving methods by Jarzynski's equality 11 VI.Numerical results for examining the passivity of thermal/ground states 12 Motivated by the progress of quantum simulation in the cold atom experiments [1][2][3][4][5] and the theoretical understanding of the thermalization of isolated quantum systems [6][7][8], the non-equilibrium dynamics of quantum systems have been studied extensively [9][10][11][12][13].A system can be driven into nonequilibrium simply by introducing a time-dependent interaction or a time-dependent coupling constant, and the system can return to equilibrium after turning off the time dependence of the coupling.For a many-body system, obtaining the exact dynamical evolution is generally difficult.Given a numerical method of evaluating the dynamical evolution, finding a way of estimating the numerical accuracy becomes a challenging arXiv:2308.13366v4[cond-mat.stat-mech]29 Apr 2024 task, especially for a system without prior knowledge of the exact dynamics.One would expect that the numerical error accumulates as the system evolves.Therefore, a reliable realtime error estimator will help compare the numerical results with the experimental ones, such as the quantum simulation.
For nonequilibrium dynamics of an isolated quantum system, one would expect that the underlying microreversibility should manifest in some way, in contrast with the macroscopic second law of thermodynamics.Indeed, equality relations exist, collectively known as fluctuation theorem [14][15][16][17][18][19][20][21][22][23], for such manifestation.The work done in a nonequilibrium process is not a state variable and will depend on the path connecting the initial and final state.The fluctuation theorem relates the average work to the free energy difference between the initial and final states.To calculate the average work or the higher moments, one can construct the corresponding characteristic function (or generating function of work statistics), which can be rewritten as the real-time correlation function [17,19,24].This characteristic function of work takes a form similar to the out-of-time-ordered correlator (OTOC) in characterizing the quantum chaos [25][26][27].In general, one needs to adopt some numerical method to calculate this generating function or some OTOC-like quantities involving nonequilibrium dynamical evaluation, which may be spoiled by the accumulation of numerical errors.On the other hand, the difference in free energy only depends on the initial and final states.It will require disproportionately less numerical effort than the one needed for the characteristic function.Exploiting the feature of such disproportionality in the numerical efforts on both sides of the fluctuation theorem, one can use it to monitor the accuracy of any numerical method for evaluating the real-time nonequilibrium dynamical evolution of a quantum system.In recent works [28], a numerical method based on the matrixproduct-state (MPS) approach [29,30] is adopted to evaluate the characteristic function for Ising spin chains and verify the fluctuation theorem, see also [31] for using the similar method for evaluating the heat statistics of open systems.This method can also be generalized to 2-dimensional spin-lattices, called the tensor-network state (TNS) [32].In this paper, we will adopt the MPS method to evaluate the characteristic function of other quantum spin chains and use the fluctuation theorem to benchmark the capability of these numerical time-evolving methods.
Another interesting aspect of work statistics is characterizing the (dynamical) quantum phase transitions (QPTs).This was firstly demonstrated in [33] for the Ising chain and later in [34] for the random matrix models under a sudden quench process, i.e, with a sudden change of the Hamiltonian, for which the work characteristic function is related to the Loschmidt echo.Moreover, the average work done on the isolated system measures the expectation value of the jump of the Hamiltonian across the quenching point with respect to the initial state.Then, a similar trick is applied to show the universal scaling behavior of the kink statistics for the dynamical quantum phase transition under the Kibble-Zurek mechanism (KZM) for the (Ising) spin chains [35][36][37][38], and for Kitave honeycomb model [39].This approach to characterizing QPTs by work statistics is extended to the (exactly solvable) spin/fermion chain models under the non-quench processes [35,40,41].
It is also known that entanglement structures like entanglement entropy can also characterize QPTs by showing a discontinuity at critical points in many cases of QPTs, e.g., see [42,43].As both the entanglement entropy and the work statistics can be encoded as some combinations of multi-point correlations, this may explain why both can characterize QPTs through some singular multi-point correlations.This perspective has been explored in [44] for the fermionic Hubbard model with random impurities.It showed that the entanglement is minimized while the work average is maximized at the critical point, but the second moment of work vanishes.In this paper, we will study this issue further for the Ising and Haldane chains.
When considering the quantum phase transition, the initial state is the ground state.After the sudden quench, it becomes a linear combination of the various excited states of the new Hamiltonian close in energy.If both initial and final Hamiltonians are non-critical and gapped, then only the gapped excited states with their energies less than the average work done on the system can be induced.However, if the final Hamiltonian happens to be at the quantum critical point, it describes a gapless system, then an infinite number of gapless excitations will be induced.We expect the work statistics to differ from those for the non-critical cases because this large number of gapless excitations will provide more dynamic paths for the quenching process for retrieving a definite amount of work.This implies that the average work excitation can display discontinuous behavior when the sudden quench connects a non-critical Hamiltonian to a critical one.In this sense, the work statistics, i.e., the average work done or the higher cumulants, can be adopted to characterize the quantum phase transitions, i.e., as an order parameter.Since the above reasoning concerns only the gapless feature of the critical systems, it can work for both the quantum phase transition of the Landau-Ginzburg type due to spontaneously-symmetry-breaking (SSB) or the topological type characterized by nonlocal order parameters.Finding the appropriate order parameter to indicate the topological phase transition is usually more difficult.In this paper, we will explore the power of the numerical MPS-based algorithm to evaluate the moments/cumulants of work statistics and demonstrate their capability and limitations in characterizing the QPTs of the quantum Ising and Haldane chains.With the help of the quantum-state RG of MPS, we can evaluate the work statistics for the chains up to 1024 sites.This will effectively suppress the finite-size effect, especially when near the quantum critical point.
Finally, we would like to study the issue of the passivity of ground states of quantum spin chains under cyclic impulse processes.It had been known that the relativistic thermal states are always passive, i.e., no work can be extracted from the system in any cyclic process [45][46][47].We first show that the passivity of thermal states is guaranteed by the fluctuation theorem or, equivalently, the second law of thermodynamics.It is also straightforward to see that the passivity of ground states under hermitian action is guaranteed by the variational principle, i.e., the ground state has the lowest energy.However, the variational principle fails to ensure the passivity of the ground states if the action for driving nonequilibrium is non-hermitian.We then adopt the MPS numerical method to examine the passivity for such cases.In all cases we checked, we found no active ground states even under non-hermitian action.Besides, the patterns of average work extraction show some interesting features.
Based on the work statistics and the fluctuation theorem, we will use reliable numerical methods to investigate the above three aspects of the quantum spin lattice models in this paper: characterizing the quantum phase transitions, benchmarking the accuracy of the numerical real-time evolution, and examining the passivity.These studies enlarge the perspective of the work done and the fluctuation theorem.In summary, in this work we will demonstrate the power of the numerical MPS method in evaluating the work statistics precisely enough to characterize the quantum phase transitions and examine the passivity of quantum states, even for the non-Hermitian actions.Moreover, we will also show how to exploit the simplicity and universality of the fluctuation theorem to benchmark the numerical real-time evolutions for generic models and numerical algorithms.
The rest of the paper is organized as follows.For completeness, in the next section, we will briefly review the basics of work statistics and the fluctuation theorem and elaborate on the theoretical frameworks of the three aspects of work statistics we will address in this paper.In section III, we will review the numerical methods for evaluating the real-time correlators based on MPS and then apply them to evaluate the generating function of work statistics and the expectation values of physical observables.In section IV, we apply the above numerical methods to Ising-like and Haldane-like spin chains.Our numerical results demonstrate that the work statistics from the sudden quench can be used to characterize the quantum phase transitions for all spin models.In section V, we will show the results of using Jarzynski's equality as the benchmark to gauge the numerical accuracy of the real-time evolving methods.In section VI, we present our results of examining the passivity of the ground states of the quantum spin chains under both hermitian and non-hermitian actions.Finally, in VII, we summarize our works and discuss the possible extensions.Some numerical consistency checks and supplements are presented in the Appendices.

A. Work statistics and fluctuation theorem
We first give a brief sketch of the basics of work statistics in the context of the fluctuation theorem, which provides an elegant framework for describing and characterizing thermal or quantum fluctuations in statistical mechanics.
Unlike some conservative quantities, work is not a state variable and will depend on the microscopic details/paths of the non-equilibrium process.In short, work is a random variable and is characterized by a distribution function defined for work statistics as follows: where p(b, t f |a) is the transition probability from the energy eigenstate |a⟩ of energy E a (0) at time t = 0 to another eigenstate |b, t f ⟩ of energy E b (t f ).This non-equilibrium process is driven by a time-dependent Hamiltonian H(t) from an initial state of density matrix ρ(0) = a p a |a⟩⟨a|, with For simplicity, in this paper, we will only consider the cases with H(t) = H 0 + λ(t)V with time-independent H 0 and V , so that [H(t), H(t ′ )] = 0 and the time-ordering symbol T in U (t f ) can be omitted.Moreover, we will in general allow for non-hermitian H(t) so that U (t f ) can be non-unitary, i.e., U † (t f ) ̸ = U −1 (t f ).The last equality implies detailed balance or micro-reversibility.Similarly, one can define the work statistics from the "reverse process" with an "initial state" of density matrix ρ(0) = m q m |m, t f ⟩⟨m, t f | as follows: where the transition probability for the reverse process is p(a|b, t f ) = | ⟨a|U † (t f )|b, t f ⟩ | 2 = p(b, t f |a) with the last equality ensured by micro-reversibility.
The fluctuation theorem is usually formulated for the canonical initial and final states, i.e., p a = e −β(Ea(0)−F (0)) and q b = e −β(E b (t f )−F (t f )) with β = 1 k B T the inverse temperature, and F (0) and F (t f ) the free energy for the initial and final states, which are defined by the partition functions Z(0) = e −βF (0) and Z(t f ) = e −βF (t f ) , respectively.Then, the fluctuation theorem taking the form of Crooks' relation [16] can be obtained by construction from Eq. ( 1)-(3), and states where ∆F = F (t f ) − F (0). Integrate this relation over W , and we can obtain Jarzynski's equality [14,17] e −βW = e −β∆F (5) where the "overline" denotes the average over work W .Using Jensen's inequality, Jarzynski's equality can yield the second law: ∆S := W − ∆F ≥ 0. The fluctuation theorem for more general quantum processes and end states can be found in [20,23,48].
From the work statistics, we can define the corresponding characteristic function as This function is the generating function of the moments of work done, e.g., For simplicity, later we will mostly adopt W = lim s→0 ∂G(−is) ∂s since G(−is) is real and G(0) = 1.Moreover, Jarzynski's equality can be expressed as follows, After some straightforward manipulation, the work characteristic function G(u) of a non-equilibrium process, which drives the initial state ρ(0) by a time-dependent Hamiltonian H(t) during the time interval [0, t f ], can be put into the following form [19], The denominator is unity for unitary U (t f ) (or hermitian H(t)).The expression of Eq. ( 10) can be recast into the realtime correlation function on the extended Schwinger-Keldysh contour [24] so that techniques of open system dynamics can be adopted.We usually omit t f and write G(u).
To be specific, in this paper, we will only consider the timedependent Hamiltonian of the following form, where we choose the initial Hamiltonian H 0 = H(0 − ) to be time-independent and hermitian; however, the driving operator V is time-independent but could be non-hermitian.If V is non-hermitian, then U (t) is non-unitary so that the denominator of Eq. ( 10) is not unity.An initial state ρ 0 , which is taken to be either the thermal state or ground state of H 0 , will be driven away from equilibrium due to the nontrivial timedependence of the coupling constant λ(t).The notations 0 − and 0 + used later denote the moments right before and after t = 0, respectively.We will consider two particular non-equilibrium processes in the following: • Sudden quench process with where Θ(t) is the Heaviside step function, and ∆λ a constant.For such a process, U (t f ) t f →0 = 1, H(0 − ) = H 0 and H(0 + ) = H 0 + ∆λV .
• Impulse process with where δ(t) is the Dirac delta function and λ a constant.Compared to a sudden quench process, the coupling constant is turned off right before and after the impulse so that H(0 − ) = H(0 + ) = H 0 .Thus, this process can be considered cyclic and will be implemented to study passivity.
In both kinds of processes, we will find that the average work done W extracted from G(u) can characterize (quantum) phase transitions.

B. Benchmark the numerical real-time evolution
Jarzynski's equality shown in Eq. ( 9) is an interesting relationship that relates a real-time correlator to the ratio of the partition functions, i.e., e −β∆F = Z(t f ) Z(0) .There is usually no analytical method to calculate the real-time correlators such as G(u) of Eq. ( 10) of a many-body system even though there might be some analytical ways of calculating the partition function.Even relying on the numerical method to evaluate the real-time correlators, the numerical errors accumulate more as the evolution continues.On the other hand, the partition function is a stationary quantity free of the accumulated error.Therefore, we can characterize the accumulated error by defining the following ratio: which is the ratio of the LHS to the RHS of Jarzynski's equality and is a function of evolution time t f .
By monitoring the deviation of R(t f ) from the unity, one can estimate the numerical error accumulation of a numerical method for real-time evolution in a first-principle way without the need to compare with the analytical or other numerical methods.For different numerical methods for real-time evolution, we can compare the deviations of their R(t f ) from the unity for different numerical methods for real-time evolution to benchmark their performances.Later, we will demonstrate benchmarking for the numerical methods adopted in this paper.

C. Work done by sudden quench and phase transition
In general, the work statistics or its characteristic function for the many-body system is difficult to evaluate for the complication of real-time dynamics.However, to characterize the phase transition, we can bypass the difficulty by just considering the work done by a sudden quench implemented by the Hamiltonian of Eq. (11) and Eq.(12).Denote the ground state of H 0 by ρ 0 ≡ |0⟩⟨0|, and the Hmailtoian after quench by H + , i.e., H + ≡ H 0 + ∆λV , then the characteristic function can be simplified as follows: which yields moments and the first few cumulants of the work done as follows, The quantity W , as shown for the sudden quench, measures the energy of the "excited state" |0⟩ with respect to H + , the quantity σ W measures the corresponding fluctuation of ∆H = ∆λV .Thus, W (or σ 2 W ) could act as a local order parameter to characterize the quantum phase transition while tuning H 0 .From Eq. ( 16) we can also calculate the higher moments W m or higher cumulants κ m for the sudden quench in terms of the expectation value of higher order V m 's.
We now elaborate on why the quantities in Eq. ( 16) and Eq. ( 17) can be used as the order parameter to characterize the phase transition.First, suppose the initial and final states are Gibbs states related by unitary evolution.In that case, the Hamiltonian can also be understood as the modular Hamiltonian, i.e., H mod = − log ρ.In this case, W can be rewritten as the relative entropy S rel [20], i.e., where ρ G −,+ are the Gibbs states before and after quench, respectively.In this case, W measures the distance between the initial and final Gibbs states.For the case we consider, the initial and final states are pure, W is no longer equivalent to the relative entropy.However, as explained below, W can still be used to distinguish pure states, especially the critical and non-critical ones.
The phase transition is usually characterized by the order parameter, which is the expectation value of some local or non-local observables.We can decompose the Hamiltonian so that Eq. ( 16) with m = 1 can be rewritten as The average work done W can be seen as the expectation value of E f + with its probability measure given by p(f ) = |⟨0|f + ⟩| 2 .We will tune some coupling in H 0 over a range covering different phases.Again, |0⟩ looks like the excited state to H + .When H 0 and H + are in the same phase, we expect On the other hand, if they are in different phases, we should not expect a sharp distribution for p(f ), and W will be quite different from E 0 + due to the incoherent average.For example, when H + is the gapless critical Hamiltonian at the quantum critical point, we expect the emergence of a dense set of degenerate gapless states denoted as {0 + c } such that the density of ground states p(0 Therefore, W can be used as an order parameter for phase transition.Similar arguments go for σ 2 W or higher moments/cumulants.Some preliminary study of p(f ) for quantum Ising chain under sudden quench by the exact diagonalization (ED) method is done in Appendix A, and the result is shown in Fig. 11 therein.

D. Passivity of a quantum state: thermal or ground states
A quantum state is called passive if there is no work extraction, i.e., W ≥ 0, under a cyclic process defined by H(t f ) = H(0).It was shown that the KMS state [49,50], i.e., thermal vacuum state for a relativistic quantum field theory, is passive [45][46][47].Even though there is no general criterion for checking the passivity of a generic quantum state, it can be argued from the fluctuation theorem as follows.For a cyclic process, we shall expect ∆F = 0 so that the fluctuation theorem becomes e −βW | cyclic = 1.By Jensen's inequality, this implies a second-law statement Therefore, the fluctuation theorem guarantees the thermal states' passivity, thus, the second law.Since the second law is established only in the thermodynamic limit, there is evidence that the thermal-like states of finite systems under some peculiar process can be active [51][52][53].
Based on the above result, it is worthy of mentioning the study in [54], which considered the passivity of the energy eigenstates of Ising chains under the action of local Hermitian operations.By adopting the eigenstate thermalization hypothesis [8,55,56], an energy eigenstate locally looks like a thermal state, e.g., the reduced density matrix of a local region is approximately a thermal one with some effective temperature.Thus, based on the passivity of the thermal state and ETH, an energy eigenstate will probably be passive under the action of local unitary operations.In [54], it was shown that all the energy eigenstates are passive under local operations for non-integrable Ising chains of finite size, while most energy eigenstates are passive for integrable Ising chains.This is consistent with the fact that ETH holds only weakly for integrable systems, see, e.g., [57][58][59].The above result could also help to explain the passivity of the ground states studied below and imply a negative effective temperature.
Intuitively, the ground states should be passive because the operation V excites the ground state to those with larger energies so that one cannot extract energy from the ground state.We can check passivity by first obtaining the average work extraction W cyclic ext ≡ −W for the cyclic process with H(t f ) = H(0) := H 0 from the characteristic function of Eq. ( 10), i.e., where Note that , so that Eq. ( 23) can be reduced to which is nothing but the energy difference between the ground state ρ 0 and the excited state Thus, in such cases, the passivity of the ground state is guaranteed if the variational principle holds.On the other hand, if U (t f ) is nonunitary, then ρ 1 ̸ = ρ 0 even though ρ 2 bears the look of excited state under non-unitary evolution.Therefore, the negativity of W cyclic ext of Eq. ( 23) can no longer be ensured by the variational principle.
We can conclude from the above discussions that the variational principle generally ensures the ground state's passivity.In this work, we will verify this by the numerical calculation and explore its possible violations by considering the cases of non-hermitian V so that the variational principle cannot be applied to guarantee the negativity of W cyclic ext of Eq. ( 23).For simplicity and numerical implementation, we will only consider the impulse cyclic processes, i.e., with λ(t) given by Eq. ( 13), for which we have This yields It turns out that we can easily implement MPS formalism for spin chain models to evaluate the right side of Eq. ( 28) numerically and then extract W impulse ext from the results.Alternatively, the average work extraction can also be obtained from Eq. ( 23)-Eq.( 25) by replacing U (t f ) with e −iλV and U † (t f ) with e iλV † .It is easy to see that W impulse ext vanishes exactly if [V, H 0 ] = 0. Otherwise, by assuming H 0 |0⟩ = 0 and taking small |λ| expansion up to O(λ 4 ), we have where we short-handed We see that the O(λ) term is absent, so that at O(λ 2 ), W impulse ext is independent of the sign of λ.Note also that the last terms in each order shown above all vanish for V † = V .Though this expression contains the terms with imaginary i factor, they are always associated with the expectation values of anit-hermitian operators, so the overall expression is a real quantity.
Besides, due to the peculiarity of the impulse process, there are some interesting features of W impulse ext .Firstly, if V = k σ a i for the spin-1/2 chain with σ a k for a = 1, 2, 3 are the Pauli matrices for the spin operator at site k.Then, i for the spin-1/2 chain, similar to but not the same as the hermitian case, e −iλV = e iλV † = k e λσ a k = k (cosh λ + σ a k sinh λ).Using this to evaluate W impulse ext and using the fact that H 0 |0⟩ = 0, we can obtain which saturates to a negative finite constant as λ → ∞.Note that the negativity is ensured by the variational principle and |⟨σ a k ⟩| ≤ 1.The peculiar features for both types of V of spin-1/2 chains are due to the highly non-adiabatic nature of the impulse process.Similar conclusions can be drawn for the spin-1 chain with some specific V composed of the sum of site-spin operators.
From the above discussion, we see that the passivity (or the negativity of W impulse ext ) is guaranteed by the variational principle in the small λ regime no matter if V is hermitian or not.Therefore, when exploring the possibility of violating the passivity of the ground states for the non-hermitian V , we will investigate the large λ regime, to which the variational principle cannot be directly applied.
Finally, we comment on the possible physical realization of a non-Hermitian operation, although our consideration in this work is somewhat academic.There are ways to realize non-Hermitian actions.One way is to drop the quantum jump term in the Lindblad master equation for the open quantum system by choosing the proper Lindblad operator, and the resultant equation describes the Heisenberg evolution of the density matrix operator by a non-Hermitian Hamiltonian, see, e.g., [60].This can be alternatively realized on quantum circuits by probability quantum computing with suitably dilated Hilbert space, see [61].Or, one can introduce the non-Hermitian interaction in the context of PT -symmetry.For the Ising chain, it is equivalent to introducing an imaginary magnetic field; see, e.g., [62].

III. NUMERICAL METHODS FOR IMAGINARY-AND REAL-TIME EVOLUTION
To calculate the work statistics, we must evaluate the associated generating function, i.e., the real-time correlator of the characteristic function, Eq. (10).Even for a quenching process, it still involves a nontrivial vacuum expectation value (VEV) as given by Eq. ( 15) or the average work done given by Eq. ( 16).This is usually a challenging task for many-body systems.To study the work statistics for various many-body quantum systems, a recent work [28] has adopted the numerical method based on matrix-product-state (MPS) [29,30] for quantum spin chains.This method is called time evolving black decimation (TEBD) [63,64] to use for the real-time evolution or evaluation of the VEV by its imaginary time evolution.Below, we will review the basic ideas of this method and mention its application for our purpose.
Moreover, Jarzynski's equality Eq. ( 9) relates a real-time correlator to the ratio of the partition functions of thermal states, which is far easier to evaluate than the former.Therefore, we can use Jarzynski's equality to gauge the accuracy of the numerical real-time evolution and serve as the benchmark for a given numerical method for real-time evolution.In the next section, we will benchmark the TEBD and the exact diagonalization (ED).

A. Matrix Product States
For the one-dimensional quantum many-body systems, the ground states could be expressed in terms of a matrix product state (MPS), which can be expressed as where number of sites of the spin chain, and A si is the on-site tensor.For example, A s k l,r are rank-three tensor for 1D MPS as shown in Fig. 1 (a).Moreover, we can also generalize this idea to represent quantum many-body operators, called matrix product operators (MPO), as shown in Fig. 1 (b).We call d s the physical dimension and χ cut the bond dimension.tTr is to sum over all indices of tensors.Depending on the situation, we should tune χ cut to capture the ground state's essential feature fully.Despite that, the χ cut generally required for a good approximation of ground states still yields a more efficient representation of the quantum ground states than the typical dimensions of d N s .The validity of MPS ansatz is due to the area-law nature of quantum entanglement entropy of ground states and shall not hold for highly excited or high-temperature states.
The MPS effectively expresses the ground states of quantum spin-chain models so that numerical methods can construct them more efficiently.We will use the time-evolving block decimation TEBD method [64] as described below to implement the imaginary time evolution operator to obtain the ground states of the two considered spin-chain models and then use them to evaluate the generating (characteristic) function for the work statistics.Furthermore, when checking the fluctuation theorem, one must prepare the initial states as Gibbs states, which are beyond primitive MPS.

B. Time-evolving block decimation
The time-evolving block decimation (TEBD) provides an efficient way to simulate time evolution.We use the TEBD method to solve for the above MPS numerically, namely, by acting the imaginary time evolution operator e −τ H on an initial state |ψ 0 ⟩ to determine the ground state and by acting the real-time evolution operator e −iHt to study the dynamics of the quantum lattice systems.
Denote the time-evolving state as |ψ t ⟩ = e −itH |ψ⟩ for an initial state |ψ⟩ with a given Hamiltonian H.For simplicity, we consider only the Hamiltonians made of arbitrary single-site and two-site terms, i.e., with nearestneighbor interactions, so that they can be decomposed as and hence the commutator H even , H odd ̸ = 0. We then break the evolution operations e −itH into a sequence of local gates using a Suzuki-Trotter expansion.For small enough δ, the Suzuki-Trotter expansion of order for e −itH can be written as where f 1 = (x, y) = xy, f 2 = (x, y) = x 1/2 yx 1/2 for first and second-order expansions, respectively.If we expand to higher-order terms, the Trotter error will decrease.The evolution operators can be expressed as a product of two-body gates using the Suzuki-Trotter expansion.The simulation of evolu-tion is achieved by updating the MPS by a sequential of alternating gate operations, i.e., alternating between e −iδHeven and e −iδH odd , as specified by Eq. (32).In this paper, we will implement this MPS-based TEBD method to evaluate the work characteristic function G(u) for extracting average work done W and also evaluate the partition functions for extracting the free energies that appear on the right-handed-side of the fluctuation theorem.

C. Thermal density matrix as MPO
As mentioned, the MPS is suitable for representing the ground states but not the highly excited states, such as hightemperature thermal states.This paper will mainly deal with the work statistics for ground states.On some occasions, we will also deal with thermal states, e.g., by checking the fluctuation theorem or the passivity of thermal states.A way to construct the thermal states based on MPS is the so-called algorithm of minimally entangled typical thermal states (METTS) [65], which adopts a Markov chain of product states to construct the thermal ensemble constrained by detailed balancing relation.
Instead of adopting METTS, a more direct way of constructing a thermal state based on MPS is to represent the density operator e −βH by an MPO.The partition function for evaluating free energy is to take a trace of this MPO.To construct such an MPO, one first prepares an initial identity MPO of bond dimension χ cut = 1 to make such an MPO.The MPO for the thermal state is then obtained by acting with the operator e −βH , which can be decomposed by TEBD and Suzuki-Trotter expansion into a sequence of a product of two-body gates, as shown in Fig. 1 (c).The partition function can be obtained by tracing out the two physical indices of the MPO e −βH as shown in Fig. 1(d).In Appendix B 1, we show the agreement on the work statistics obtained from METTS and MPO representation of thermal states.In the rest of the paper, we will adopt MPO to represent thermal states to proceed with the numerical calculations.

IV. NUMERICAL RESULTS FOR CHARACTERIZING THE QUANTUM PHASE TRANSITIONS
In this section, we will present our numerical results to demonstrate that the average work done by the sudden quench can be the order parameter for the quantum phase transitions.The quantum systems we consider include quantum spin-1/2 and spin-1 chains.The quantum phase transition considered for the spin-1/2 chain is the Landau-Ginzburg type due to spontaneous symmetry breaking.On the other hand, the ones for the spin-1 chain model are the topological phase transitions, which can only be characterized by non-local order parameters.
The average work done by sudden quench is given succinctly by Eq. ( 16) or extracted from the characteristic function provided by Eq. (15).In either case, we can implement the TEBD method based on MPS to evaluate it.Below, we will show the numerical results case by case.

A. Quantum spin-1/2 chain in magnetic field
We first consider the anisotropic Heisenberg XY spin-1/2 chain in a transverse magnetic field h z and a longitudinal field h x .Its Hamiltonian is given by where S x,y,z k are the site spin-1/2 operators, k labels the number of sites.When γ = 1, it reduces to the quantum Ising chain, and we will focus on the parameter range: 0 ≤ γ ≤ 1 and h x = 0.The ground state of this model can be exactly solved and exhibits three phases of the Landau-Ginzburg-Wilson paradigm: the oscillatory phase (O), the ferromagnetic phase (F ), and the paramagnetic phase (P ) [66].The ground state stays in the O phase for small h z until h z = 1 2 1 − γ 2 and then changes to the F phase as a crossover transition.Further increasing h z up to h z = 0.5, the ground state will change to P phase, which is a second-order quantum phase transition in the thermodynamic limit with ⟨S x ⟩ as the order parameter, i.e., ⟨S x ⟩ is nonzero in F (or O if γ = 0) phase but is zero in P phase.Below, we will consider the phase transition between the ferromagnetic and paramagnetic phases.W , (g) κ3 and (h) κ 4 .The plots are obtained by the TEBD method with the parameters χcut = 16 and δ = 0.01, and W n are obtained using Eq. ( 16) by the MPO method.As expected, the higher moments or cumulants yield sharper change near the critical point so that the work statistics can indicate the quantum phase transitions in the Ising-like spin chains.To go beyond sudden quench and for comparison, we also extract W and W 2 (black dashed lines) from the characteristic function Eq. 10 for the non-sudden quench process with t f = 0.1.
Here, we apply the MPS-based TEBD method to calculate the characteristic function G(u) of Eq. (10) for this spin chain model with γ = 1 and h x = 0 at zero temperature but varying the transverse field h z .We aim to check if the work done W can be used as an order parameter for the quantum phase transition, i.e., from F (or O if γ = 0) phase to P phase.For simplicity, we consider the sudden quench process implemented by the Hamiltonian of Eq. ( 11) and Eq. ( 12) with a chosen jump of the coupling constant, i.e., ∆λ = 0.01.The results are shown in Fig. 2 and 3.
The Fig. 2 shows the characteristic function G(−is) of a n = 1024 Ising chain for s ∈ [−0.2, 0.2] with γ = 1 and h x = 0 at various values of h z in the initial Hamiltonian H 0 .For simplicity, we pick five values of s but vary h z continuously with its value indicated by the color bar attached aside.In particular, the five points belonging to h z = 0.5 are denoted by the crosses and joined into a line by Lagrange interpolation.We see that the slope of this line is negative, which implies that the average work done W = lim s→0 ∂G(−is) ∂s is negative.For other values of h z in the initial Hamiltonian H 0 , we can see that the slope of the line connected by the five points of the same h z starts with a positive value for h z = 0, then gradually turns negative as h z increases and finally converges to a constant negative value.
The above results and further related calculations can be translated into phase diagrams as shown in Fig. 3 for the quantum phase transition from the F (or O) phase to the P phase.These phase diagrams are characterized by the entanglement entropy S, the average work done, W , and its associated higher moments and cumulants, under a sudden quench process with V = k S z k and ∆λ = 0.01.We see that S has a sharp peak near h z = 0.5, which implies a stronger correlation near the critical point at h z = 0.5.
Similarly, the work statistics denoted by W and the higher moments and cumulants also show the shaper changes near h z = 0.5.As expected, the higher moments and cumulants show sharper changes at the critical point.This implies that the work statistics can indicate the phase transitions.Interestingly, the cumulants of even orders show a level-off for the F (or O) phase with h z ≤ 0.5, and similar behavior for σ 2

W
has been observed earlier obtained by the exact solution of Ising chains [67].This constancy could be associated with some underlying symmetry.However, the quantities ⟨V n ⟩ with V = k S z k , related to the cumulants of work, are different from the order parameter ⟨ k S x k ⟩ for the transition from F (or O) phase to the P phase.
In addition, in Appendix B 2, we give a consistency check to show in Fig. 13 that the above numerical results of W and W 2 obtained from either Eq. ( 7) or Eq. ( 16) agree.

B. Anisotropic quantum spin-1 chain
To show that our proposal can also work for the quantum phase transitions of the non-Landau-Ginzburg-Wilson type, we consider a particular model of such type, the anisotropic quantum XXZ spin-1 chain described by the following Hamiltonian:  W , κ ℓ=3,4 and S, the subfigure captions and the numerical implementing methods are the same as the one in Fig. 3.The quantum critical points are indicated by the sharper behavior of S, which are at (i) (Jz, D) = (1.5, 0.48), (ii) (Jz, D) = (1.5, 1.34) and (iii) (Jz, D) = (1.0,1.0).Note that (i) is the transition point from the ordered phase to the topological phase, but (ii) and (iii) are the purely topological ones without any order parameter and should be more difficult to detect.For (i), W m=1,••• ,4 show very mild behavior, but cumulants show shaper changes as the order of cumulants goes higher.However, for (ii) and (iii), only σ 2 W and κ4 show mild crossover (see the enlargement insets in the subfigures (f) and (h)), and all the other moments/cumulants in this figure show no obvious changes.Overall, we can conclude that the work statistics can barely capture the purely topological phase transitions.The plots are obtained by the TEBD method with the parameters χcut = 64 and δ = 0.01.
where S x,y,z k are the site spin-1 operators, and the parameter D denotes the uni-axial anisotropy.When J z = 1, it reduces to the so-called Haldane chain.The ground-state phase diagram of this model consists of six different phases [68,69].Here, we focus on phase transitions among three phases: the large-D phase, the N'eel phase, and the Haldane phase characterized by nonzero string-order parameters.At large D, the model is in a trivial insulator phase.Namely, the ideal large-D phase is |000 • • • ⟩.On the other hand, the ideal N'eel phase is The Haldane phase is one of the symmetry-protected topological (SPT) phases.The ground states of such phases have nontrivial patterns of quantum entanglement.They cannot continuously connect to trivial product states without closing the gap or breaking the protecting symmetry.Thus, the SPT phases preserve the global symmetry of the Hamiltonian so that some spontaneous symmetry-breaking local order parameters cannot characterize the transition from the SPTnontrivial phase to the SPT-trivial one.This is in contrast to the cases of the Landau-Ginzburg type.Instead, the topological phases, such as the Haldane phase, are characterized by fractionalized edge excitations, which some nonlocal order parameters, such as the expectation value of some string-like operator, can measure.
As before, we can characterize the quantum phase transitions of this model by entanglement entropy S, the average work done W and its higher moments/cumulants by a sudden quench with V = k (S z k ) 2 and ∆λ = 0.01, and by the associated higher moments and cumulants.The results are shown in Fig. 4, which are again obtained numerically using the MPS-based TEBD method.We see that the entanglement entropy S shows the peaks around three quantum critical points, which are at (i) (J z , D) = (1.5, 0.48), (ii) (J z , D) = (1.5, 1.34) and (iii) (J z , D) = (1.0,1.0).The critical point (i) is the transition point from the ordered Neel phase to the topological Haldane phase and is partly with local order parameters.On the other hand, the critical points (ii) and (iii) are from the Haldane phase to the large D phase, and no local parameter exists to characterize it.Therefore, these transitions are purely topological types and should be more difficult to detect.Indeed, this is what we see from Fig. 4. For (i), W m=1,••• ,4 show very mild behavior, but cumulants show shaper changes as the order of cumulants goes higher.However, for (ii) and (iii), only σ 2 W shows mild crossover, and all the other moments/cumulants in this figure show no obvious changes.It is strange why the higher cumulants do not show sharper behavior for critical points (ii) and (iii).This could be due to the need for a higher bond dimension to capture the behavior of higher cumulants.Overall, we can conclude that the work statistics can barely capture the purely topological phase transitions.Since work statistics involve all multipoint correlations, and could be used as some non-local order parameter.Thus, we expect the higher moments and cumulants of the work statistics to display shaper behavior around the topological phase transitions, although it is quite computationally costly to obtain higher moments and cumulants due to the need for large bond dimensions.Also, in Appendix B 2, we give a consistency check to show in Fig. 14 that the above numerical results of W obtained from either Eq. ( 7) or Eq. ( 16) agree.However, the W 2 obtained from either Eq. ( 7) or Eq. ( 16) can't match well.The larger bound dimension χ cut might improve the physical properties near critical points.

V. BENCHMARKING THE REAL-TIME EVOLVING METHODS BY JARZYNSKI'S EQUALITY
We now adopt the Jarzynski equality of Eq. ( 9) dictated by the real-time correlator expression of the work characteristic function G(iβ; t f )| β=1 , i.e., Eq. ( 10), to benchmark the numerical accuracy of the real-time evolution of two quantum spin chain models considered above.The numerical error mainly comes from evaluating G(iβ; t f ) and accumulates as t f grows.For the small-size cases, we can use either the Exact Diagonalization (ED) method or the MPS-based TEBD method.For large-size cases, we can only use the TEBD method.As before, the non-equilibrium process is implemented by the Hamiltonian of Eq. (11), and here, we adopt the linear profile of coupling constant, i.e., λ(t) = t.The results for the small-size cases are shown in Fig. 5 for (a) n = 10 Ising chain of fixed h x,z and with V = i S z i and (b) n = 6 Haldane chain of fixed D = 2 and with V = i (S z i ) 2 .Fig. 5 0.0 0.2 0.4 0.6 0.8 1.0 Benchmarking the ED and TEBD methods by the ratio k (S z k ) 2 .In subfigure (c), we show the dependence of the resultant R(t f ) of n = 6 Haldane chain on the bond dimension χcut of MPS.As expected, larger χcut yields better results as expected.The results from ED are indicated by squares, and from TEBD, they are indicated by circles.Besides, we also show the ratio Z(0) |ED (triangles in subfigures (a) and (b)) to demonstrate the accuracy of TEBD in evaluating the partition functions.In all the above, the coupling constant has the linear profile, i.e., λ(t) = t.If there is no numerical error, R = 1, and the deviation benchmarks the real-time numerical errors.For ED, we see that R remains equal to one.On the other hand, for MPS, we see that the real-time numerical errors start to escalate around t f = 0.7.However, the ratio of partition functions remains accurate.
shows that as t f increases, the curve of G(iβ; t f )| β=1 of the Ising chain increases; on the other hand, the curve in the Haldane chain decreases.
Next, we use Jarzynski equality to characterize the numerical error of real-time evolution and benchmark the corresponding numerical methods.This can be quantified by the ratio R(t f ) = G(iβ;t f ) Z(t f )/Z(0) | β=1 defined in Eq. ( 14), which is the ratio between either side of Jarzynski equality.If there is no numerical error, the Jarzynski equality is obeyed so that R = 1.Otherwise, the |1 − R| size can characterize the numerical error.
We can evaluate the partition function ratio Benchmarking the numerical accuracy with Jarzynski equality for n = 100 Ising chain with the same H0, V , χcut and δ as in Fig. 6(a).The results for G(iβ; t f )| β=1 , R(t f ) and are denoted by blue circles, red crosses, and black triangles, respectively.Due to the large size, ED is unavailable, and the results are obtained by the TEBD method.We see that the real-time numerical errors start to escalate around t f = 0.13, which is far shorter than the n = 10 case as expected.However, the ratio of partition functions remains accurate as in the n = 10 case.
ED and TEBD methods.As expected, both methods can be accurate for such non-dynamical quantity.Combined with the previous results for G(iβ; t f )| β=1 , we can evaluate R. The results for small-size quantum chains are shown in Fig. 6 with subfigure (a) for n = 10 quantum Ising chain and subfigure (b) for n = 6 Haldane chain.We see that R remains one for the ED method.However, for the MPS-based method, the numerical errors escalate around some critical moment, i.e., (t f ) c = 0.7 as shown in Fig. 6 (a) and (b).Moreover, we also evaluate and show the ratio of the ratios of the partition functions obtained by ED and MPS-based methods, i.e., | ED (black triangles in Fig. 6 (a) and (b)).In Fig. 6 (c), we show the dependence of the resultant R(t f ) of n = 6 Haldane chain on MPS's bond dimension χ cut .The larger χ cut yields better results as expected.In this case, χ cut ≥ 40 is needed to arrive at an accurate R = 1 result for t f ≤ 1.
When the quantum spin chain size becomes large, the ED method is unavailable, and we can only use the MPS-based method to evaluate R. For comparison, we evaluate the benchmarking factor R for n = 100 Ising chain with the same H 0 , V , χ cut and δ as in the n = 10 case of Fig. 6(a).The results of G(iβ; t f )| β=1 , R(t f ) and are shown in Fig. 7.
As we can see, G(iβ; t f )| β=1 grows exponentially with t f , which should cause large numerical errors.Indeed, the realtime numerical errors characterized by |1−R| start to escalate around the critical moment (t f ) c = 0.13.This is far shorter than (t f ) c = 0.7 of the n = 10 case, as we expect that the numerical errors accumulate more quickly for systems of larger sizes.
. Demonstration of passivity of thermal states in n = 10 quantum Ising chain under impulse process implemented by the Hamiltonian of Eq. ( 11) and Eq. ( 13) with H0 of Eq. ( 33) with hx = hz = 0 and with V = k S x k , k S y k or k S z k .The resultant average work done per site w impulse for β = 1 are obtained by ED and TEBD.Its positiveness for all considered V 's and λ ∈ [0, 0.5] demonstrates the passivity of thermal states as pointed out in Eq. ( 22).

VI. NUMERICAL RESULTS FOR EXAMINING THE PASSIVITY OF THERMAL/GROUND STATES
We now can consider the passivity of a given initial state undergoing a (cyclic) impulse process.Two natural candidate categories are thermal or ground states.The passivity of (relativistic) thermal states was shown to be ensured [45][46][47].For the thermal states of the non-relativistic systems, such as the quantum spin chains, as argued in Eq. ( 22), the passivity should be guaranteed by the fluctuation theorem, i.e., the second law of thermodynamics.Here, we demonstrate this is the case for the n = 10 Ising chain.The result is shown in Fig. 8, from which we see that the average work done per site w impluse = −w impulse ext for the thermal state at β = 1 is always positive for various V and λ.
We now turn to the issue of passivity for the ground states under the cyclic impulse processes.As shown in section II D, the passivity of ground states is guaranteed to be passive under the Hermitian action, i.e., V † = V by the variational principle.This is because the average work extraction from a cyclic process can be understood as the energy difference between the ground state and the excited state driven by V .It is always negative as the variational principle guarantees.On the other hand, there is no such simple interpretation for the average work extraction if the action V is non-Hermitian, i.e., one cannot apply the variational principle to check the passivity of ground states.For the latter cases, it is then interesting to examine the passivity by evaluating average work extraction numerically.
As a startup, we numerically examine the passivity of ground states for various Hermitian or non-Hermitian actions but small λ, i.e., λ = 0.1.In Fig. 9, We plot the average work done per site w impluse = −w impulse ext as the function of h z , which labels the ground states of the Ising-like chain under various Hermitian and non-Hermitian actions.A similar result for the Haldane-like chains is also given in Fig. 16 of Appendix C. To obtain these results, we first evaluate G(−is) by the MPS-TEBD method, and then from it extract w impluse = −w impulse ext from G(−is).The ground states considered are all passive under the various Hermitian and non-Hermitian actions described in Fig. 9.This is a consistency check for our numerical method with the discussion in section II D for examining the passivity of ground states of Ising-like chains.
By exploiting the numerical method's power in checking the ground states' passivity, we now consider the cases with quite an extensive range of λ ∈ [−10, 10] for both Hermitian and non-Hermitian actions.The results for the Ising-like chains are shown in Fig. 10, and the similar results for the Haldane-like chains and the more general non-Hermitian actions in the Ising-like chains are shown in Fig. 17 and Fig. 15 of Appendix C, respectively.Since the ones shown in the Appendix C bear similar features as in Fig. 10, we will not discuss them in the main text.The top row of Fig. 10 shows the average work done per site on the ground states labeled by the values of h z ∈ [0, 1] for the Hermitian actions V = k S x,y,z k , and the bottom row shows the results for the non-Hermitian actions V = i k S x,y,z k .For the Hermitian cases, we see a quite interesting feature that W impulse is periodic with respect to the coupling λ with the period equal to 2π or π.This periodic behavior has been discussed in section II D due to the peculiar feature of the impulse process.We will not see such behavior for the cyclic process with a finite time duration.For the non-Hermitian case, we first see that all the ground states considered are passive as in the Hermitian cases.However, W impulse is no longer periodic with respect to λ nut levels off as λ becomes large.In both cases, the average work done is bounded no matter how large the coupling λ is.Again, this is due to the zero time duration of the impulse process, so the average work done is limited to such a tiny time duration.However, it is unclear why the passivity remains intact for non-Hermitian action, especially without the guarantee by the variational principle.This issue deserves future study to explore.

VII. CONCLUSIONS
In this paper, we apply the numerical method for realtime evolution, such as Exact Diagonalization (ED) and MPSbased TEBD, to quantum spin lattice models to study the work statistics.With the power of MPS formalism, we can evaluate the chains up to 1024 sites, which can effectively suppress the finite size effect.We focus on three aspects: (i) study the behaviors of the moments and cumulants of the work statistics near the quantum phase transitions and examine their capability to indicate quantum critical points; Our results up to the fourth cumulant show that the work statistics can detect the quantum phase transitions charactered with local order parameters, but just barely for the topological phase transitions; (ii) we propose to adopt Jarzynski equality as the benchmark for the accuracy of the numerical real-time evolution methods; (iii) we examine the passivity of thermal states and ground states of quantum spin chains under some cyclic impulse processes.Our numerical results show that all the ground states are passive under both Hermitian and non-Hermitian actions considered in this work.Although the variational principle ensures the passivity of ground states under Hermitian actions, it is not the case for non-Hermitian actions.It is interesting to explore the reason for the passivity of ground states under non-Hermitian action seen in this paper, and also explore the possibility of active ground states by more general actions.Once the active ground states exist, we may adopt them to implement the quantum engine naturally to extract quantum work in the cyclic processes.
The quantum spin lattice with nearest neighbor interactions can be the natural system for performing quantum simulation and can serve as accurate tests for work statistics such as fluctuation theorem.Due to the statistical nature of quantum work, it has remained quite mysterious since its proposal decades ago.With our demonstration of numerical studies for the realistic many-body systems, one can explore more different perspectives of work statistics.For example, when using pure states as the initial states for the work statistics, it may need more subtle treatment than the two-point measurements to preserve the quantum coherence and explore its role in the fluctuation theorem.We may hope to adopt the real-time evolution method, such as the MPS-based one, to investigate their role in some specific quantum tasks.21) are plotted as a function of the eigen-energy E f of the Hamitain H+ = H0 + ∆λV for the n = 10 quantum Ising chains with hz = 0.4, 0.5, 1.0 under a sudden quench process with V = k S z k and ∆λ = 0.1.We can see that the p(f )'s for hz = 0.5 (black crosses) and hz = 0.4 (red circles) have a more focused energy spectrum than the one for hz = 1.0 (green pluses).Note that for hz = 0.5 (hz = 0.4), the initial (final) Hamiltonian is in the critical phase because ∆λ = 0.1.Thus, the more focused energy spectrum for these two cases implies a tendency for the effect of ground-state degeneracy.
12. Numerical agreement of the results for the Jarzynski equality.We evaluate work characteristic function G(iβ; t f )| β=1 (crosses) and the ratio of partition functions (circles) for n = 20 quantum Ising chain by TEBD and MPO of χcut = 32 and δ = 0.01 to represent the thermal states.The result shows Jarzynski equality holds and justifies using MPO for constructing the thermal states.After comparison, we also find that our results numerically agree with the ones from the ensemble of 10000 METTS, as shown in Fig. 4 in Ref. [28].
This can justify our usage of MPO representation of thermal states for examining the passivity of thermal states, as shown in Fig. 8.
2. Equivalence between average works done obtained from Eq. ( 10) or from Eq. ( 16) There are two ways to obtain the average work done by a sudden quench process.The first one is to extract it from the work characteristic function G(u) of Eq. ( 10) by using Eq. ( 7).The second one calculates it by evaluating the vacuum expectation value (vev) of the Hamiltonian offset using Eq. ( 16).We check the agreement of the results from both ways for the quantum Ising and Haldane chains of n = 10 by ED or n = 1024 by MPS-based method, as shown in Fig. 13 and Fig. 14, respectively.
In Fig. 13 for the quantum Ising model with various h z in the initial Hamiltonian H 0 under sudden quench process with V = k S z k and ∆λ = 0.01, we show that the numerical agreement of the numerical results of W m for m = 1, 2, 3, 4 obtained from the characteristic function G(−is) (red crosses) or from the vev of V m (blue circles).This serves as a consistency check for the TEBD method.
In Fig. 14 for n = 1024 Haldane (XXZ) chain with various coupling D in the initial Hamiltonian H 0 under sudden quench process with V = k (S z k ) 2 and ∆λ = 0.01, we show that results for (a) W , (b) W 2 , (c) σ 2 W and (d) entanglement entropy S, which are obtained from the characteristic function G(−is) (black dashed line for χ cut = 32 or solid line for χ cut = 64) or from the vev of V m (red crosses for χ cut = 64 or blue circles χ cut = 64).We can see that the results from both ways match well for (a) and (c) but not for (b).This mismatch can be due to the need of large χ cut to compensate for the errors of performing numerical derivatives on G(−is).This is seen from the improvement of W 2 in (b) and S in (d) when enlarging χ cut twice large.shows the plot of W impulse as a function of λ to examine the passivity of the ground states under either the Hermitian or non-Hermitian actions of a cyclic impulse process.The key features are similar to what we have observed in Fig. 9 and Fig. 10 of the main text for the Ising chains.In particular, the ground states are all passive for all the cases considered.Moreover, the features shown in Fig. 17 are kind of the mixtures of the top and bottom rows in Fig. 10 as the V 's considered here are the mixtures of those twos.We see that the results from both ways match well for (a) and (c) but not for (b).We thus evaluate entanglement entropy S in (d) to show the relevance of increasing the bond dimension χcut in improving the accuracy, which could also be the reason for the mismatch in (b).
statistics and fluctuation theorem 3 B. Benchmark the numerical real-time evolution 4 C. Work done by sudden quench and phase transition 4 D. Passivity of a quantum state: thermal or ground states 5 III.Numerical methods for imaginary-and real-time evolution 6 A. Matrix Product States 7 B. Time-evolving block decimation 7 C. Thermal density matrix as

FIG. 1 .
FIG. 1.(a) Symbol for matrix product state (MPS) structure.(b) Symbol for matrix product operator (MPO) structure.(c) MPO representation of the operator e −βH for a quantum system with nearestneighbor interaction via TEBD algorithm.(d) MPO representation of approximated e −βH .

FIG. 2 .
FIG. 2. Work characteristic function G(−is) under a sudden quench process implemented by Hamiltonian H = H0 + ∆λ V with ∆λ = 0.01 with H0 of a n = 1024 quantum Ising chain, i.e., (33) of γ = 1, hx = 0 and hz ∈ [0, 2], and V = k S z k .The results are obtained by the MPS-based TEBD method via the 1-st order Trotter expansions for δ = 0.01 with bond dimension χcut = 16.We show the results only for five values of s but can connect them into a line by Lagrange interpolation.Each color represents one particular choice of hz, with their values indicated by the color sidebar.The Lagrange interpolating lines at critical points are shown explicitly for hz = 0.5.The slope of an interpolating line corresponds to the average work done W , which is mostly negative.

with tf = 0 6 FIG. 3 .
FIG. 3. Moments and cumulants of work done by sudden quench with V = k S z k and ∆λ = 0.01 on quantum Heisenberg XY Ising chain of size n = 1024 as a function of the transverse field hz at zero temperature for γ = 0.5 (red crosses), 0.8 (blue circles) and 1.0 (green squares).For comparison, the entanglement entropy S is also shown.They are listed as follows: (a) W , (b) W 2 , (c) W 3 , (d) W 4 , (e) S, (f) σ 2W , (g) κ3 and (h) κ 4 .The plots are obtained by the TEBD method with the parameters χcut = 16 and δ = 0.01, and W n are obtained using Eq.(16) by the MPO method.As expected, the higher moments or cumulants yield sharper change near the critical point so that the work statistics can indicate the quantum phase transitions in the Ising-like spin chains.To go beyond sudden quench and for comparison, we also extract W and W 2 (black dashed lines) from the characteristic function Eq. 10 for the non-sudden quench process with t f = 0.1.

Jz = 1 6 FIG. 4 .
FIG.4.Moments and cumulants of work done by sudden quench with V = k (S z k )2 and ∆λ = 0.01 on quantum Heisenberg XXZ Haldane chain of size n = 1024 as a function of the coupling parameter D at zero temperature for JZ = 1.0 (red crosses) and 1.5 (blue circles).For comparison, the entanglement entropy S is also shown.The listing order for W m=1,••• ,4 , σ 2 W , κ ℓ=3,4 and S, the subfigure captions and the numerical implementing methods are the same as the one in Fig.3.The quantum critical points are indicated by the sharper behavior of S, which are at (i) (Jz, D) = (1.5, 0.48), (ii) (Jz, D) = (1.5, 1.34) and (iii) (Jz, D) = (1.0,1.0).Note that (i) is the transition point from the ordered phase to the topological phase, but (ii) and (iii) are the purely topological ones without any order parameter and should be more difficult to detect.For (i), W m=1,••• ,4 show very mild behavior, but cumulants show shaper changes as the order of cumulants goes higher.However, for (ii) and (iii), only σ 2 W and κ4 show mild crossover (see the enlargement insets in the subfigures (f) and (h)), and all the other moments/cumulants in this figure show no obvious changes.Overall, we can conclude that the work statistics can barely capture the purely topological phase transitions.The plots are obtained by the TEBD method with the parameters χcut = 64 and δ = 0.01.

FIG. 5 .
FIG. 5. Work characteristic function G(iβ; t f )| β=1 implemented by Hamiltonian H = H0 + λ(t)V with λ(t) = t, for (a) n = 10 Ising chain, and (b) n = 6 Haldane chain, with their H0 and V shown on the corresponding subfigures.We obtain the results by using the ED method (circles) and the MPS-based TEBD method (crosses) with the parameters χcut = 40 and δ = 0.01.In some cases, both results from ED and MPS agree.Note that G(iβ; t f )| β=1 increases with time for (a) but decreases for (b).

from
Jarzynski equality for subfigure (a) n = 10 Ising chain with hx = hz = 1 and V = k S z k ; and subfigure (b) n = 6 Haldane chain with D = 2 and V =

5 SFIG. 9 .
FIG.9.Plot of w impulse to show the passivity of ground states of n = 1024 quantum Ising chain with hx = 0 but varying hz ∈ [0, 1] under impulse process implemented by Hamiltonian of Eq. (11) and Eq.(13) with λ = 0.1 and Hermitian actions V = k S x k (red circles ), k S y k (black cross), k S z k (green cross), and non-Hermitian actions V = k iS x k (blue diamonds), k iS y k (yellow square), k iS z k (pink triangles).Again, we see that W impulse can also indicate quantum phase transition by mild crossover behaviors when compared with the phase diagram of the entanglement entropy of initial state (gray circle).

FIG. 10 .
FIG.10.Plot of w impulse to show the passivity of ground states and its pattern as a function of the coupling λ of n = 1024 quantum Ising chain.We fix hx = 0 but vary hz ∈ [0, 1] to obtain various ground states, and also vary λ ∈ [−10, 10] under impulse process implemented by Hamiltonian of Eq. (11) and Eq.(13) with the Hermitian actions (a)V = k S x k , (b) V = k S y k , (c) V = k S z k inthe top row, and the non-Hermitian actions (d) V = i k S x k , (e) V = i k S y k and (f) V = i k S z k in the bottom row.

FIG. 11 .
FIG.11.The ED results of weight p(f ) from Eq. (21) are plotted as a function of the eigen-energy E f of the Hamitain H+ = H0 + ∆λV for the n = 10 quantum Ising chains with hz = 0.4, 0.5, 1.0 under a sudden quench process with V = k S z k and ∆λ = 0.1.We can see that the p(f )'s for hz = 0.5 (black crosses) and hz = 0.4 (red circles) have a more focused energy spectrum than the one for hz = 1.0 (green pluses).Note that for hz = 0.5 (hz = 0.4), the initial (final) Hamiltonian is in the critical phase because ∆λ = 0.1.Thus, the more focused energy spectrum for these two cases implies a tendency for the effect of ground-state degeneracy.
Appendix C: More results for examining the passivity of ground states in Ising-like and Haldane-like chainsIn this Appendix, we present more results of examining the passivity of the ground states.The first one is shown in Fig.15for the plot of W impulse of Ising-like chains as a function of λ ∈ [−10, 10] for the more general non-Hermitian actions.It can be seen as the non-Hermitian counterpart of Fig.10.Then, we present the Haldane-chain counterparts of Fig.9and Fig.10, respectively, in Fig.16and Fig.17 .
Fig.16shows the plot of W impulse as a function of coupling parameter D labeling the ground states of the Haldane-like chains under various Hermitian actions for small value of λ = 0.1.Fig.17

4 FIG. 13 .FIG. 14 .
FIG.13.Agreement of W m of m = 1, 2, 3, 4 obtained from either characteristic function G(−is) of Eq. (10) (red crosses) or from the vev of V m (blue circles) for the transverse quantum Ising chains under sudden quench process with V = k S z k and ∆λ = 0.01.In (a) and (b), the result for n = 10 chains is obtained by ED method, and for (a) it also shows the result from the vev of Hamiltonian offset by Eq. (16) (green squares).In (c),(d),(e), and (f), the results for n = 1024 chains are calculated by the MPS-based method.

FIG. 15 .
FIG.15.Plot of w impulse to show the passivity of ground states and its pattern as a function of the coupling λ of n = 1024 quantum Ising chain.We fix hx = 0 but vary hz ∈ [0, 1] to obtain various ground states, and also vary λ ∈ [−10, 10] under impulse process implemented by Hamiltonian of Eq. (11) and Eq.(13) with the non-Hermitian actions (a) V = k (S x k + iS y k ), (b) V = k (S y k + S y k ),and (c) V = k (S x k + S z k ).

FIG. 16 .
FIG.16.Plot of w impulse to show the passivity of ground states of n = 1024 Haldane-like chain with Jz = 1.5 but varying D ∈ [0, 2] under impluse process of λ = 0.1 and Hermitian ac-tions V = k S x k (black triangles), k S y k (yellow right-arrows), k S z k (red circles), k (S x k ) 2 (green diamonds), k (S y k ) 2 (orange left-arrows) and k (S z k ) 2 (blue squares).The W impluse plots show that the ground states are always passive for all considered V 's.The degeneracy of W impluse under swapping of S x i and S y i , as shown here, can be understood as the Z2 symmetry of H0 and SO(3) spin algebra under the transformation: (S x k , S y k , S z k ) → (−S y k , S x k , S z k ).Besides, W impluse always vanishes for V = k S z k because of [H0, k S z k ] = 0.

2 FIG. 17 .
FIG. 17. Plot of w impulse to show the passivity of ground states and its pattern as a function of the coupling λ of n = 1024 Haldane-like chain with Jz = 1.5.We vary the coupling D ∈ [0, 2] to obtain various ground states, and also vary ∈ [−10, 10] under impulse process implemented by Hamiltonian of Eq. (11) and Eq.(13) for the Hermitian actions (a)V = k S x k , (b) V = k (S x k ) 2 , (c) V = k (S z k )2 in the top row, and non-Hermitian actions (d)V = i k S x k , (e) V = k (S x k + iS y k ), and (f) V = i k (S x k ) 2 in the bottom row.