Closed hierarchy of Heisenberg equations in integrable models with Onsager algebra

Dynamics of a quantum system can be described by coupled Heisenberg equations. In a generic many-body system these equations form an exponentially large hierarchy that is intractable without approximations. In contrast, in an integrable system a small subset of operators can be closed with respect to commutation with the Hamiltonian. As a result, the Heisenberg equations for these operators can form a smaller closed system amenable to an analytical treatment. We demonstrate that this indeed happens in a class of integrable models where the Hamiltonian is an element of the Onsager algebra. We explicitly solve the system of Heisenberg equations for operators from this algebra. Two specific models are considered as examples: the transverse field Ising model and the superintegrable chiral 3-state Potts model. Copyright O. Lychkovskiy. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 22-03-2021 Accepted 17-05-2021 Published 01-06-2021 Check for updates doi:10.21468/SciPostPhys.10.6.124


Introduction
Describing an out-of-equilibrium dynamics of a quantum many-body system is, in general, a formidable task. One can expect that this task is considerably simplified for integrable models. However, the integrability turns out to be somewhat deceptive when it comes to the dynamics. A range of methods for addressing the dynamics and resulting non-thermal steady states is under active development, including the quench action approach [1,2], generalized Gibbs ensemble [3], generalized hydrodynamics [4,5] and advanced techniques for summing formfactor expansions [6][7][8]. Nevertheless, each instance of an analytical calculation in this field is a certain tour de force.
In the present paper we address the out-of-equilibrium many-body dynamics via a system of coupled Heisenberg equations. In a generic many-body model, such system of equations consists of exponentially many equations involving progressively nonlocal operators, analogously to the closely related Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy of equations on reduced density matrices [9][10][11][12][13]. However, one can expect that in an integrable system a subset of operators can turn out to be closed with respect to the commutation with the Hamiltonian, and, as a result, the Heisenberg equations for these operators decouple from the rest of the equations. This indeed happens for a quadratic fermionic (or bosonic) Hamiltonian, where for any fixed integer k a set of operators spanned by products of k fermionic (or bosonic) creation and annihilation operators remains closed with respect to commutation with the Hamiltonian. This property of quadratic Hamiltonians can also be used to describe the dynamics spin-1/2 models that can be mapped to systems of noninteracting fermions. For example, in Ref. [14] an open-end transverse field Ising model was considered, and a set of Heisenberg equations for operators linear in terms of fermionic operators (k = 1) was derived and solved.
Here we study models with the Hamiltonian of the form where a 0 , a 1 are real numbers and A 0 , A 1 are two operators that satisfy the Dolan-Grady conditions [15] This pair of operators generates the Onsager algebra of operators [16] (as briefly reviewed below in Sec. 2.1). This class of models includes the transverse-field Ising model [17][18][19][20] and superintegrable chiral Potts models [21]. The latter are not reduced to noninteracting fermions. We derive and solve a system of Heisenberg equations for a set of operators G n , A n spanning the Onsager algebra. This is done in the next section. Then, in Sec. 3, we apply the general solution to describe the out-of-equilibrium dynamics of the transverse field Ising model for a translation-invariant product initial state polarized in an arbitrary direction. In Sec. (4) we address the out-of-equilibrium dynamics of the 3-state Potts model. We conclude by the summary and outlook in Sec. 5.

Heisenberg equations and their solution 2.1 Onsager algebra
The Onsager algebra [16] is a Lie algebra spanned by operators G n , A n that are recursively generated starting from A 0 , A 1 according to Commutation plays the role of the bilinear product of the algebra, with the structure explicitly given by with G −n = −G n . The Dolan-Grady conditions (2) are necessary and sufficient for the set of operators A n , G n generated by the recursion (3) to satisfy eq. (4) [15,22,23]. Originally, the Onsager algebra emerged in studies of the classical Ising model [16]. An apparently first explicit construction of an Onsager algebra representation for a quantum model -namely, for the X Y spin chain -was presented in [24]. In an important paper by von Gehlen and Rittenberg [21] it was shown that for an arbitrary positive integer n a nearest-neighbour spin-n/2 Hamiltonian of the form (1) can be constructed, with A 0 and A 1 satisfying the Dolan-Grady conditions (2) and thus generating a representation of the Onsager algebra. The spin-1/2 case is the transverse-field Ising model. The higher spin models are known as n-state superintegrable chiral Potts models or Z n -invariant clock models. This discovery triggered a considerable lasting interest in such models [25][26][27][28] and, more generally, in applications of the Onsager algebra and its generalizations to quantum integrability [29][30][31][32][33].

Solving Heisenberg equations
We will work in the Heisenberg representation where a time-dependent expectation value of an observable O is given by with ρ 0 being the initial state of the system and O t -the Heisenberg operator related to the Schrödinger operator O as The Heisenberg operator satisfies the Heisenberg equation of motion, with the initial condition O 0 = O.
We are going to solve Heisenberg equations for Heisenberg operators G n t , A n t . These equations are straightforwardly obtained from eqs. (1), (4): where n = 0, ±1, ±2, . . . We exclude A n t from eq. (8) and obtain a system of the second-order equations: This system can be conveniently written in the compact matrix notations as where G t is a vector composed of G n t and T is a tridiagonal Toeplitz matrix with matrix elements T nn = −16(a 2 0 + a 2 1 ) and T n n±1 = −16 a 0 a 1 . To proceed further, we truncate the matrix T to a finite (N − 1) × (N − 1) matrix, keeping the same notation for the truncated matrix. In fact, this truncation emerges naturally for finite-dimensional representations of Onsager algebra [26], where A n , G n are operators on the Hilbert space of a spin chain with N spins, as discussed in more detail in the next section. Anyway, the dependence on the system size goes away for local observables in the thermodynamic limit of N → ∞.
We then diagonalize the truncated matrix by a unitary transformation U, where the matrix elements of U read and E is a diagonal matrix with diagonal entries By a standard argument, the diagonalization (11) implies the following solution of eq. (9): with where the limit N → ∞ is already taken. A n t is then found from the second line of eq. (8) and reads where Eqs. (14)- (17) are the main general results of the present paper. They express timedependent Heisenberg operators G n t , A n t through the Schrodinger operators G n , A n .

Out-of-equilibrium dynamics of the Ising model
In the present section we apply our general results to the transverse field Ising model. This model belongs to the simplest class of noninteracting integrable spin chains since it can be mapped to noninteracting fermions [17][18][19][20]. Early studies of the transverse field Ising model mostly addressed equilibrium instantaneous and time-dependent correlation functions [14,20,[34][35][36][37][38][39][40] (this topic still attracts a deal of attention [8,41,42]). The studies of the outof-equilibrium dynamics, while initiated in notable early papers [34,43], came into focus much later [8,[44][45][46][47][48][49][50], with most efforts directed to the dynamics after a quantum quench, as reviewed e.g. in [51]. Here we describe the out-of-equilibrium dynamics for a previously unstudied initial condition specified below. We consider the transverse field Ising chain with N spin 1/2. The Hamiltonian is given by eq. (1), where the finite-dimensional representation of the Onsager algebra reads (see e.g. [52]) Here σ x, y,z j are Pauli matrices, n = 1, 2, . . . , (N −1), all operators G n , A n (as well as the Hamiltonian) are assumed translation invariant with σ α N +n ≡ σ α n , and n−1 m=1 σ z j+m is void (i.e. replaced by 1) for n = 1.
Note that, thanks to the equality G N = 0, the first (N − 1) equations in eq. (9) decouple from the subsequent equations, which rigorously justifies the truncation of the matrix T in eq. (10). This also implies that our results, though presented here in the thermodynamic limit, are also valid, upon a straightforward modification, for any finite N .
Let us remark that the "string" operators (18) recurrently show up in the studies of the Ising model [20] and related models, the topics varying from integrable Floquet dynamics [53] to a variational ansatz for nonintegrable spin chains [54] to adiabatic gauge potentials [55].
We consider the dynamics of the Ising model prepared in a product initial state where p = (p x , p y , p z ) is the polarization vector inside the Bloch sphere, p 2 ≤ 1. While for the initial polarisation along the z-axis the problem has been studied earlier [34], for a generic polarisation the dynamics has remained unexplored. We wish to calculate time-dependent expectation values 〈G n 〉 t , 〈A n 〉 t . To this end we plug the solution (14) and . We plug these values into eq. (14) and explicitly take the sum over m, which is essentially a geometric series. This way we obtain in the thermodynamic limit where In Fig. 1 we plot the evolution of We have verified that in a special case of p z = ±1, p x = p y = 0 our result for 〈σ z j 〉 t coincides with that of Niemeijer [34]. Further, very recently the dynamics starting from the state with p x = 1, p y = p z = 0 has been studied in detail in ref. [56]. An analytical formula for 〈σ x j σ x j+1 〉 t provided in [56] agrees with our result. Let us emphasise that the results (21) have been obtained without resorting to the fermionic picture of the Ising model. Undoubtedly, the same results could be obtained in the fermionic language, since the string operators (18) are just quadratic operators in terms of fermions [24].  (31). Plotted is the difference ∆ l between two approximations to 〈S z j 〉 t derived from the truncated eq. (31), with truncation orders equal to l max = l and l max = 4, respectively. The parameters of the Hamiltonian are a 0 = 1.2 and a 1 = −1. One can see that truncation error is small already for l max = 1 and rapidly vanishes with growing l max .
However, our method is more straightforward, both conceptually and technically. Furthermore, it avoids known difficulties of the fermionic approach related to the boundary term dependent on the overall parity j σ z j [39], as well as to dealing with initial states not amenable to the generalized Wick's theorem [51].
We also note that for equilibrium time-dependent correlation functions a set of nonlinear equations has been derived [57,58] and solved numerically [59]. This method relies on the Wick's theorem and, therefore, can not be generalized to the dynamics starting from non-Gaussian initial states.

Out-of-equilibrium dynamics of the 3-state Potts model
Here we consider the superintegrable chiral 3-state Potts model [21]. It is a model of N spins-1 with the Hamiltonian (1), where Here ω = e 2πi/3 and τ j , σ j are operators acting on the j'th spin with the following properties: Analogously to the previous section, the system is assumed to be translation invariant with τ j+N ≡ τ j , σ j+N ≡ σ j . This model, in contrast to the Ising spin-1/2 chain considered above, does not map to noninteracting fermions.
We choose the following matrix representation of these operators (see e.g. [28]): With this choice, the z-projection of a spin, S z ≡ diag(1, 0, −1), can be expressed as S z = ((1 − ω * ) −1 τ + h.c.), and therefore A 0 is, up to the factor 4/3, a total polarization along the z-axis: Importantly, τ j and τ † j do not change the z-polarization of the j'th spin, while σ j , σ j τ j , σ † j τ j and their conjugates change the polarization by one. We will refer to the latter type of operators as shifting operators.
In contrast to the Ising model, an explicit closed form of G n , A n is not known [28]. This means that the exact expressions (14), (16) for Heisenberg operators G n t , A n t can not be immediately converted to exact expressions for expectation values 〈G n 〉 t , 〈A n 〉 t , as in the case of the Ising model. However, it turns out that sums in eqs. (14), (16) converge so rapidly that it suffices to calculate a few first G n , A n to obtain an excellent approximation to the exact result. Below we demonstrate this by calculating 〈S z j 〉 t for a simple initial pure state ρ 0 = |Ψ 0 〉〈Ψ 0 | polarized along the z-axis, where | ↑〉 is the spin-up state, S z | ↑〉 = | ↑〉. We employ computer algebra to calculate the first few G n , A n starting from A 0 , A 1 according to the recursive relations (4). The resulting expressions have a considerably more complicated appearance than those in the Ising case. For example, The number of essentially different terms necessary to represent higher G n , A n rapidly grows: it equals 9 for A 2 and A −1 , 24 for G 2 , 43 for A 3 and A −2 , 100 for G 3 , 181 for A 4 and A −3 , 424 for G 4 etc. Inspecting these expressions, we conjecture the following properties of G n , A n .
1. Each G n , A n is a sum of strings, where a string of length m is a tensor product of m operators, each chosen from the set acting on m consecutive sites.
2. Each string of length n ≥ 2 has shifting operators at both its ends.
3. Each G n consists of strings of lengths not less than 2.
4. Each A n with odd n consists of strings of lengths not less than 2.
5. Each A n with even n contains strings of length 1 consisting of operators τ j , τ † j .
We leave the proof (or refutation) of these conjectures for further work. In actual calculations we use only first few G n , A n where these properties have been verified explicitly. From these properties it follows that 〈G n 〉 0 = 0 for all n (in fact, this can be easily derived from the first of the recurrence relations (3)) and 〈A n 〉 0 = 0 for odd n. Further, for even n the only type of string that contributes to 〈A n 〉 0 is the string of length 1 constructed from τ j or τ † j . Finally, from the second line of eq. (4) with m = 0, one obtains 〈A n 〉 0 = 〈A −n 〉 0 . As a result, we obtain from eq. (16) a simple expression with C 1(2l−1) t given by eq. (17). Note that a 0 enters this formula through eqs. (13), (17). In contrast to the Ising case, an explicit general formula for 〈A 2l 〉 0 is not known. For this reason we have to truncate the sum in eq. (31) at a finite l = l max . Fortunately, the convergence is very rapid. The first few 〈A 2l 〉 0 calculating iteratively from eq. We plug these values into eq.(31) and approximately compute 〈S z j 〉 t , see Fig. 2 (a). We compare four truncations in eq. (31) with l max = 1, 2, 3, 4 and observe a very rapid convergence, as illustrated in Fig. 2 (b).
We conclude this section by emphasizing that a successive application of our method is conditioned to the understanding the Onsager algebra representation for a specific model under consideration. A progress in such understanding in highly desirable.

Summary and outlook
To summarize, we have solved coupled Heisenberg equations for a set of operators comprising an Onsager algebra in a system with the Hamiltonian of the form (1), which itself belongs to this algebra. The solution is given by eqs. (14)- (17). It allows one to obtain the time evolution of the corresponding observables as soon as the expectation values of these observables in the initial state are known. We have considered two specific realization of the Hamiltonian (1). The first one describes the transverse-field Ising model, where we have calculated time dependence of an infinite set of observables for a translation-invariant product initial state (19) with an arbitrary polarisation, see eqs. (21) and Fig. 1. As a second example, we have studied the superintegrable chiral 3-state Potts model. There we have obtained an approximate but highly accurate description of time evolution of the transverse polarization, see eq. (31) and Fig. 2.
Our approach can be extended to algebras different from the Onsager one. In particular, a set of ∼ 4N 2 strictly local operators entering the sums in eq. (18) is closed with respect to commutation. Therefore, one can attempt to derive and solve Heisenberg equations for individual strictly local operators from this set. This will allow one to study the dynamics for non-translation-invariant initial states and/or transverse field Ising Hamiltonians (in particular, inhomogeneous quantum quench setups [60]) in a site-resolved manner.
An interesting open question is whether the method presented here can be extended to a broader range of integrable model. Most integrable models are not known to possess a simple algebraic structure analogous to the Onsager algebra (see, however, a recent ref. [61] where a hidden Onsager algebra has been conjectured for the integrable XXZ spin-1/2 chain at the root-of-unity anisotropies). The absence of such structure prevents a straightforward generalization of the method.
Note, however, that one does not necessarily need an algebra of operators to apply the method. Indeed, the requirement of closeness of the set of operators with respect to the commutation is excessive: One actually needs a more weak property of being closed with respect to commutation with the Hamiltonian. Further work will show whether the method can be applied beyond the integrable models with the Onsager algebra.