Group-theoretical approach to the calculation of quantum work distribution

Usually the calculation of work distributions in an arbitrary nonequilibrium process in a quantum system, especially in a quantum many-body system is extremely cumbersome. For all quantum systems described by quadratic Hamiltonians, we invent a universal method for solving the work distribution of quantum systems in an arbitrary driving process by utilizing the group-representation theory. This method enables us to efficiently calculate work distributions where previous methods fail. In some specific models, such as the time-dependent harmonic oscillator, the dynamical Casimir effect, and the transverse XY model, the exact and analytical solutions of work distributions in an arbitrary nonequilibrium process are obtained. Our work initiates the study of quantum stochastic thermodynamics based on group-representation theory.

Usually the calculation of work distributions in an arbitrary nonequilibrium process in a quantum system, especially in a quantum many-body system is extremely cumbersome. For all quantum systems described by quadratic Hamiltonians, we invent a universal method for solving the work distribution of quantum systems in an arbitrary driving process by utilizing the group-representation theory. This method enables us to efficiently calculate work distributions where previous methods fail. In some specific models, such as the time-dependent harmonic oscillator, the dynamical Casimir effect, and the transverse XY model, the exact and analytical solutions of work distributions in an arbitrary nonequilibrium process are obtained. Our work initiates the study of quantum stochastic thermodynamics based on group-representation theory.

PACS numbers:
Introduction.-In the past two decades, a great breakthrough in the field of nonequilibrium thermodynamics is the discovery of fluctuation theorems [1] and the establishment of stochastic thermodynamics [2,3]. It extends the usual definition of work, heat and entropy production in a thermodynamic process from ensembleaveraged quantities to trajectory-dependent quantities. Based on these extensions, the second law is sharpened and rewritten into equalities (fluctuation theorems) for arbitrary nonequilibrium processes. For an isolated quantum system, a proper definition of the trajectory work is defined through the so-called two-point measurement [4,5], which preserves the non-negativity of the probabilities of trajectories and Jarzynski equality at the same time [6,7]. As is known, the work distribution of a system provides a great deal of information about a nonequilibrium process [8][9][10], which is an analogue to the partition function encoding essential information about an equilibrium state. However, very few analytical results of work distributions of quantum systems have been obtained so far. In literature, the few exceptions are associated with either specific models (e.g., a forced harmonic oscillator [11][12][13], a harmonic oscillator with a time-dependent frequency [14,15], quenched Luttinger liquid at zero temperature [8], a diving scalar field [16] or a transverse Ising model at zero temperature [15]) or very special protocols (e.g., the sudden quench protocol [17][18][19][20][21][22][23] or the quantum adiabatic protocol [24,25]). Hence, to develop a universal method to efficiently calculate the work distribution for various quantum systems in an arbitrary nonequilibrium process is one of the most challenging problems in this field.
In this letter, we invent a universal method which is based on the Lie-group theory to calculate the work distribution for quantum systems described by quadratic Hamiltonians [26,27]. Such Hamiltonian is ubiquitous in physics. Examples include the harmonic oscillator, some lattice models, phonons and Bloch electrons in solids, identical particles in an external potential, some models in quantum field theory [16] and many other physical models within the mean-field approximation, e.g., superconductors and superfluids [27]. The results of work distributions can be solved exactly and efficiently for an arbitrary work protocol. And in many cases, the analytical solution can be obtained explicitly, which not only has pedagogical value, but also brings important insights to the design of optimal protocols for thermal machines [28]. Moreover, our method reduces the cost of computational resource from the exponential growth [29] to the polynomial growth with the particle number. Since the work distribution provides a great deal of information about thermodynamics in a nonequilibrium process, our method will hopefully significantly advance the development of quantum stochastic thermodynamics, especially for quantum many-body systems.
Two-point measurement scheme.-In this scheme, the trajectory work for an initial canonical ensemble is determined by the projective measurements over the Hamiltonians before and after the work protocol. Suppose the system evolves from time t = t 0 to time t = t 1 with the initial stateρ and the time-dependent Hamiltonian H(t). For simplicity, letρ commute withĤ(t 0 ). We measure the energy of the system at the initial and final time over the HamiltoniansĤ 0 =Ĥ(t 0 ) andĤ 1 =Ĥ(t 1 ) respectively. Also, we assume the states of the system are projected into |E 0 m and |E 1 n after the corresponding measurements. Here, {|E 0 m } and {|E 1 n } are the instantaneous orthonormal bases, i.e.,Ĥ 0 |E 0 m = E 0 m |E 0 m and H 1 |E 1 n = E 1 n |E 1 n . Then the trajectory work for the above realization is given by w nm = E 1 n − E 0 m . After repeating the above procedure with the same canonical ensemble and the same driving protocol numerous times, the work distribution is obtained as follows whereÛ (t 1 , t 0 ) = T exp[−i t1 t0 dsĤ(s)] is the unitary evolution operator (we set = 1 in this letter) and T is the time-ordered operator. For theoretical analysis, it is a better choice to consider the Fourier transform of the work distribution which is known as the characteristic function of work (CFW) G(v) [30] where the index H inĤ H 1 indicates that the operator is in Heisenberg picture, i.e.,Ô H (t) =Û † (t, t 0 )Ô(t)Û (t, t 0 ). It is worth mentioning that G(v) encodes the same information as P (w).
Calculating method of the CFW.-As long as the Hamiltonian is a quadratic polynomial of bosonic or fermionic operators and the initial state is a thermal equilibrium state, the CFW is solvable. In the following, we elucidate our method for arbitrary work protocols.
For later convenience, we introduce some notations first. In a Hilbert space, . . , 2n. denote the annihilation and creation operators for bosons (fermions). These operators can be considered as the elements of two column vec- Now, without losing generality, we first consider a time-dependent bosonic system whose Hamiltonian can be written aŝ where A(t) is a 2n × 2n matrix, B(t) is a 2n-dimensional column vector and C(t) is a scalar. Here, A(t), B(t) and C(t) are work parameters that are controlled by an external agent. From Eq. (2), we see that the calculation of the CFW can be decomposed into two steps. The first step is the calculation ofĤ H 1 , and the second step is to calculate the trace of the product of several exponential operators. Based on these facts, the CFW for systems described by the Hamiltonian (Eq. (3)) can be obtained according to the following procedure: . . ,â †H n (t)) T denote α in Heisenberg picture. Then α H (t) satisfies the following equation Because the Hamiltonian is quadratic, the differential equation (Eq. (4)) is linear. As a result, the solution of α H (t) can be written in the following linear form where D(t) is the general solution and E(t) is the particular solution of the linear ordinary differential equation (Eq. (4)) [31]. Substituting Eq. (5) into Eq. (4), we ob-tainḊ where the overhead dot denotes the time derivative, τ B is a constant matrix given by and I (I) denotes the 2n × 2n (n × n) identity matrix.
, we obtain the expression of the Hamiltonian in Heisenberg picture. Please note that the HamiltonianĤ H 1 is still a quadratic polynomial of α.
(2) Calculate the trace of the product of several exponential operators in Eq. (2) by utilizing the Lie-group representation technique which will be elucidated in next section.
Similarly, we can also consider a time-dependent fermionic system whose Hamiltonian can be written aŝ We would like to emphasize that when there is no linear part in Eq. (7), i.e., B(t) ≡ 0, the calculating method for quadratic fermionic systems is almost the same as that for quadratic bosonic systems mentioned above. However, when B(t) ≡ 0, the above method fails, because the Bogoliubov transformation does not apply in this situation [27,[32][33][34]. For this situation, we can employ a new method which is described as follows to bypass this difficulty.
We introduce an artificial Hamiltonianˆ H(t) by adding a pair of identical fermionic 'ghost' operatorsb 0 and T is linear again. Now, the artificial Hamiltonian has no linear part. The above method can be applied, and we obtain the CFW G(v) corresponding toˆ H(t) by similarly following the procedure for bosons. Finally, the CFW corresponding toĤ(t) is obtained by Lie-group representation technique.-A key point in the calculation of the trace (Eq. (2)) is to calculate the product of several exponential operators and this can be significantly simplified by utilizing the Lie-group representation technique. A number of different representations of Lie algebras can be constructed from the bosonic (fermionic) annihilation and creation operators [35][36][37]. For bosons, we have (1 denotes the identity operator) (ii) The 2n-dimensional symplectic algebra sp(2n) (iii) The semidirect product of symplectic algebra and Heisenberg algebra sp(2n) s h(n) And for fermions, we have where δ ij denotes the Kronecker delta function.
(v) The 2n + 1-dimensional special orthogonal algebra so(2n The above algebras are all over the complex number field C. We would like to emphasize that all quadratic Hamiltonians can be constructed by the above Lie algebras [35]. Hence, our study has exhausted [38] all situations of quantum systems described by quadratic Hamiltonians. The exponential map maps the Lie algebras to the corresponding Lie groups, the Heisenberg group H(n), the symplectic group Sp(2n), the semidirect product of the symplectic group and Heisenberg group Sp(2n) s H(n) and the spin groups Spin(2n) and Spin(2n+1) [35]. Thus, the product (group multiplication) of several exponential quadratic operators with or without a linear part is equal to another exponential quadratic operator with or without a linear part [37]. The trace of the latter can be obtained [39][40][41] by using the coherent-state representation technique for bosons and fermions. In the following, we give the trace formula corresponding to the product of several exponential quadratic operators without a linear part.
(1) Quadratic bosonic operators without a linear part: where m denotes the number of the exponential operators.
(2) Quadratic fermionic operators without a linear part: where are the 2n-dimensional representation matrices of Sp(2n) (Spin(2n)), acting on the vector space spanned by the basis The symmetric (antisymmetric) condition of S i (R i ) is due to the commutation (anticommutation) relation of bosons (fermions) [37]. Also, we require the operators in the trace are trace class operators and τ B S i (τ F R i ) are invertible and diagonalizable [26].
For quadratic bosonic operators with a linear part, we can also use the coherent-state representation and Weyl ordering rule [44] to obtain the trace formula: (3) Quadratic bosonic operators with a linear part: Sp(2n) s H(n) where l 1 and l 2 are 2n-dimensional column vectors. Here L reads where For quadratic fermionic operators with a linear part, it is a bit tricky. By adding a pair of ghost operators (see supplemental material) [34] and using Eq. (10), we obtain the trace formula (4) Quadratic fermionic operators with a linear part: Spin(2n + 1) ⊂ Spin(2n + 2) where R i is determined by the following equation We would like to emphasize that the sign of the square root of a complex number in Eqs. (9)(10)(11)14) is determined by two properties of the CFW: G(0) = 1 and G(v) is a smooth function of v. Moreover, some trace formulas corresponding to the subgroups of the Lie groups mentioned above are specials cases of our results. For example, the Levitov's formula [42,43] which corresponds to unitary groups [38] can be reproduced from Eqs. (9,10).
Example.-In order to illustrate the effectiveness of our method, in the following, we solve the CFW of a 1D generic time-dependent harmonic oscillator which belongs to the group Sp(2n) s H(n) .
Let us consider a generic non-homogeneous secondorder linear ordinary differential equation This equation describes a forced harmonic oscillator with a time-dependent frequency ω(t) and mass m(t). In addition, it is subject to an external driving force m(t)c(t).
From the equation of motion (Eq. (16)), one can derive the Hamiltonian Then the parameter a(t) in Eq. (16) can be expressed as For simplicity, we set c(t 0 ) = 0.
To quantize the system, we employ the canonical commutation relation [x,p] = i. Then the solutions ofx H (t) andp H (t) read with y 1 (t 0 ) = 1,ẏ 1 (t 0 ) = 0, y 2 (t 0 ) = 0,ẏ 2 (t 0 ) = ω(t 0 ), f (t 0 ) = 0, where y 1 (t) and y 2 (t) are the general solutions and is the particular solution of Eq. (16). In the derivation, we have used the Abel's identity [31] where m 0 = m(t 0 ) and ω 0 = ω(t 0 ). Also we define the Fock states asâ †â |n = n |n , n = 0, 1, . . .. Thus, we can rewrite the Hamiltonian into a quadratic polynomial of operators,â andâ † , and obtain the CFW by the calculating method mentioned above. As for the initial state, we would like to calculate the CFW when the system is initially prepared in a microcanonical state or a canonical equilibrium state. For this purpose, we introduce an artificial density operatorρ(s) = e lnsâ †â = n s n |n n| where s is a parameter. Then, the nth Fock state (the microcanonical state)ρ n = |n n| and the canonical equilibrium stateρ th = e −βĤ0 /Z 0 , Z 0 = Tr(e −βĤ0 ) can be generated fromρ(s) bŷ where β is the inverse temperature of the initial thermal state. Using Eq. (11), we obtain the generating function of the CFW G(v, s) when the initial state is chosen to bê ρ(s)

Group-theoretical approach Previous results
Harmonic oscillator with a time-dependent frequency Generic time-dependent (when m(t) = m0, c(t) = 0) [14] Boson harmonic oscillator: Sp(2n) s H(n) Forced harmonic oscillator (arbitrary protocol: m(t), ω(t), c(t))(Eq. (23)) (when m(t) = m0, ω(t) = ω0) [11] Quenched Luttinger Liquid at zero temperature [8] A driving quantum scalar field [16] Transverse XY model: Spin(2n) Transverse XY model Fermion (arbitrary protocol) (Eq. (S28)) (sudden quench protocol) [17][18][19] Anisotropic XY model in a magnetic field parallel to X-or Y-axis: Spin(2n + 1) [34,46] where m 1 = m(t 1 ), ω 1 = ω(t 1 ) and c 1 = c(t 1 ). Thus, the CFW for the microcanonical initial state G n (v) and for the canonical initial state G th (v) can be generated from G(v, s) by It is straightforward to check that the above expressions of the CFW, G n (v) and G th (v), satisfy the following conditions: (1) G n (0) = G th (0) = 1; (2) G th (iβ) = Z 1 /Z 0 , where Z 1 = Tr(e −βĤ1 ). The two conditions indicate the normalization condition of the CFW and the validity of Jarzynski equality. We would like to emphasize that our results (Eq. (23)) include the results of the forced harmonic oscillator [11] and the result of the harmonic oscillator with a time-dependent frequency [14] as special cases. That is, for a forced harmonic oscillator, m(t) = m 0 , ω(t) = ω 0 , we reproduce Eqs. (30,33) in Ref. [11] (S + is denoted as |z| 2 in that reference); Ror a harmonic oscillator with a time-dependent frequency, m(t) = m 0 , c(t) = 0, we reproduce Eq. (17) in Ref. [14] (R + is denoted as Q in that reference). As for other examples, the CFW associated with the dynamical Casimir effect is obtained in Ref. [45] which belongs to the group Sp(2n). And the CFW of a 1D Transverse XY model is given in the supplemental material which belongs to the group Spin(2n). In order to demonstrate the powerfulness of our method, we compare the results obtained from the Lie-group representation technique with previous results in  [37,39], we invent a universal method based on Lie-group theory to calculate the trace of products of several exponential quadratic operators. Our method is valid as long as the Hamiltonian can be written in the form of the quadratic polynomial of bosonic or fermionic operators which is ubiquitous in physics especially in quantum many-body systems. Our trace formulas exhausted all possible Lie groups for quadratic Hamiltonians. Hence, it is the maximum extension of Levitov's formula [42,43]. By utilizing these trace formulas, we develop a systematic and general procedure of calculating the work distributions in many quantum systems under arbitrary work protocols in a much more efficient way. In comparison with previous results which are restricted to either several specific models or very special protocols, we can obtain the analytical expression of the CFW by a unified method. In addition, our method reduces the cost of the computational resource from the exponential growth (the dimension of Hilbert space) to the polynomial growth (the dimension of the representation space of Lie group) with the particle number. As examples, we calculate the CFW of a generic time-dependent harmonic oscillator and the CFW of a time-dependent transverse XY model (see supplemental material) to illustrate the powerfulness of our method. We also would like to emphasize that this method has potential applications in many other fields. Such as quantum optics and statistical physics. Examples include the partition functions, Out-of-time-ordered correlation functions [47][48][49], Loschmidt echos [50,51], Full-Counting Statistics [52,53] and other correlation functions. Studies about applications of our method in these problems will be given in forthcoming papers. Let H denote the 2 n -dimensional Hilbert space of a fermionic system with the basis vectors |{n i } ≡ ⊗ i |n i , n i = 0, 1 and the ladder operatorsb j ,b † j (i, j = 1, 2, . . . , n). We introduce a 2 n+1 -dimensional Hilbert space H by adding a pair of fermionic ghost operatorsb 0 ,b † 0 with the basis vectors |0 , |1 . These operators obey the anticommutation relation where {Û ,V } =ÛV +VÛ . Then, the Hilbert space H can be divided into two 2 n -dimensional subspaces H = H + ⊕ H − with the basis vectors |{n i } + ⊕ |{n i } − , where Correspondingly, the orthogonal projective operators on the two subspaces areP ± = 1 2 [1 ± (b 0 +b † 0 )], where1 denotes the identity operator. According to Ref. [1], we construct an isomorphic mapping between the Hilbert space H and the subspace H + in the following way:  where we have used the propertyP 2 where x, y, z are the Cartesian coordinates,σ x,y,z l are Pauli operators, J is the coupling strength, γ is the anisotropic parameter and Γ denotes the strength of the external magnetic field in the z direction. Also, we assume the Born-von Kármán cyclic condition and the lattice points are labeled with l = 1, 2, . . . , N, N + l ≡ l. We assume N to be a even number for simplicity.
Also, the above expression of the CFW G(v) satisfies the normalization condition of the CFW, G(0) = 1 and Jarzynski equality, G(iβ) = Z 1 /Z 0 . If the number of the lattice is very large, i.e., N → ∞, the energy spectral of the two subspaces will be the same [2,8]. Hence we have .

(S32)
Also, we obtain the first two moments of work We would like to emphasize that from our main results Eq. (S28), one can derive the CFW in the sudden quench limit and quantum adiabatic limit: 1. Sudden quench limit X k (t 1 ) = 1 0 0 1 .
Hence, we have Q k = cos[φ k (t 1 )−φ k (t 0 )]. Then, the expression k∈K + g + k (v) is exactly Eq. (11) in Ref. [7] which is the CFW of the 1D Transverse XY model in the sudden quench limit without considering the negative-parity Hamiltonian.

Quantum adiabatic limit (N is finite)
Hence, we have Q k = 1, which is an indication of quantum adiabaticity (see Eq. (S30)).