Quantum kinetic Ising models

We introduce a quantum generalization of classical kinetic Ising models, described by a certain class of quantum many body master equations. Similarly to kinetic Ising models with detailed balance that are equivalent to certain Hamiltonian systems, our models reduce to a set of Hamiltonian systems determining the dynamics of the elements of the many body density matrix. The ground states of these Hamiltonians are well described by matrix product, or pair entangled projected states. We discuss critical properties of such Hamiltonians, as well as entanglement properties of their low energy states.


Introduction
The paper that we have prepared for the special issue of the New Journal of Physics on "Quantum Information and Many-Body Theory" has two apparently independent motivations. First, it is motivated by the recent interest in many body quantum master equations and design of open systems for quantum state engineering and quantum simulations. Second, it concerns entanglement properties of quantum many body states corresponding to classical kinetic models with detailed balance and their generalizations to quantum master equations.

Quantum master equations
The quantum master equation (QME) is the basic theoretical tool to describe the evolution of open systems that undergo Markovian dynamics [1,2]. While the initial studies of QMEs dealt with many body applications, over the last decades most research has focused on systems of few degrees of freedom, mainly because of the unquestionable complexity of the many body QMEs. In the recent fifteen years, however, there have been two waves of renewed interest in many body QMEs, related to the unprecedented progress of experimental control and engineering of ultracold atomic and molecular systems. On one hand, there has been a wealth of interest in the studies of various kind of cooling processes using QMEs. Zoller, Gardiner, and collaborators [3] studied in a series of papers the growth of a Bose-Einstein condensate in a trapped bosonic gas. One of us (M.L.), together with I. Cirac, P. Zoller, Y. Castin, and others, developed the theory of laser cooling of Bose gases [4], in particular in the so called festina lente limit [5]. Similar ideas were applied to the processes of laser cooling of Fermi gases [6] and sympathetic cooling [7,8,9].
On the other hand, several authors proposed to make use of the capabilities of modern quantum optics and atomic physics experimental methods to design and realize experimentally QMEs that allow to engineer interesting quantum many body pure states [10,11,12]. These pure states range from simple Bose Einstein condensates to the stabilizer states, MPS-PEPS states, or states with topological order. The first experiment realizing these ideas was conducted in the group of G. Rempe [13], who has been able to prepare a 1D bosonic Tonks-Girardeau [14] gas employing threebody losses. So far most of these new proposals concern ultracold atoms, molecules, in particular polar molecules, or ions, but recently they start to involve Rydberg atoms, that are particularly suitable for design and realization of 3-body, 4-body etc. interactions thanks to the Rydberg blockade mechanism. Weimer et al. [15], proposed to use many-body quantum gates stroboscopically, employing the Rydberg blockade effect to engineer the topologically ordered ground state of the famous Kitaev toric code, the color code [16], and even to realize a quantum simulator of the U (1) lattice gauge theory.
Motivated by these developments we propose in this paper a new class of many body QMEs generalizing classical kinetic models (in particular kinetic Ising models (KIMs); for a review see [17]). Our QMEs have as stationary states thermal Boltzmann-Gibbs states of the underlying classical model. Note that these thermal states might correspond to very complex quantum states, for instance if the underlying classical model concerns kinetics of commuting stabilizer operators for the cluster states [18], as we discuss in Sec. 5. The diagonal elements of the density matrix in our models undergo dynamics equivalent to that of the underlying classical model. The off-diagonal elements of the density matrix exhibit complex evolution, but fortunately can be grouped into independently evolving blocks. Similarly to kinetic Ising models with detailed balance, which are equivalent to certain Hamiltonian systems, our models reduce to a set of Hamiltonian systems determining the dynamics of the elements of the many body density matrix. In 1D we identify classes of these models that are exactly soluble. In general, the ground states of these Hamiltonians are well described by the matrix product states (MPS), or pair entangled projected states (PEPS) [19,20].
It should be noticed that there exist previous efforts in the above direction. For instance, in Ref. [21] (see also Ref. [17]) a general master equation for kinetic Ising model with an environment was derived. Under some assumptions, the diagonal elements of these formulations also reproduce Glauber's kinetic model, although no full solution of the equation was attempted.

Entanglement in many body systems
The studies of the role of entanglement in many body systems were definitely initiated by the seminal Ref. [22], but they go back to the early works on area laws in quantum systems [23] (for a review see [24,20,25]. Let us remind the readers that grounds states of (non-critical) many body systems with local Hamiltonians exhibit the area law. This means that if we divide the whole system into a subsystem A and the rest B where the size |A| of A is large, but much smaller that that of |B|, and calculate the von Neumann entropy of the reduced density matrix of A, the latter will scale as the size of the boundary ("area") of A, S A ∝ |∂A|. The area law expresses the fact that away from criticality, correlations -and in particular quantum correlations responsible for entanglement-decay on short length scales.
The situation is different at criticality, although so far only 1D systems are fully understood. In one dimension, the area law at criticality may get logarithmic corrections so that the entropy of the block of size L scales as S L ∝ c log(L), where the constant c can be related to the charge of the conformal field theory describing the corresponding critical behavior. In higher dimensions the situation is much less clear, and no universal laws at criticality are known [20]. Very recently, Masanes [26] proved under quite general conditions that for the ground states of systems with local Hamiltonians, the entropy of a block A is bounded from above by |A| log A.
There is a class of quantum states that always fulfill the area law in any dimension, despite the fact that they exhibit criticality for certain values of parameters [20]. These states are related to thermal states of classical Hamiltonians, such as Ising, or Potts models. In fact, any set of local quantum mechanical commuting operators can be used to build such models. For example, one can take stabilizer operators for cluster states [18], or star and plaquette operators of the Kitaev model [16]. In the following we will mainly focus on Ising models, leaving the discussion of more complicated cases to further publications. For Ising models, the states in question have the form ‡ where σ = (σ 1 , . . . , σ N ) denotes a configuration of N Ising spins (|σ stands for a vector representing σ in the corresponding Hilbert space), H(σ) is the corresponding classical Hamiltonian, and Z N = Tr exp[−βH(σ)] is the partition function. As pointed out in Refs. [19,20], and references therein, these states have the following properties: • They are associated to a family of classical kinetic Ising models (KIMs) that describe the approach to the classical thermal equilibrium state exp[−βH(σ)]/Z N and obey detailed balance.
• The KIMs in question can be transformed by an appropriate ansatz into an equivalent problem of Hamiltonian dynamics in imaginary time (i.e. describing decay in time), with a certain quantum HamiltonianĤ parametrically depending on β. • The HamiltonianĤ is non-negative and has one eigenvector given by the expression (1) corresponding to the eigenvalue zero. Away from criticalityĤ is gapped, i.e. excited states have strictly positive energies. As β → β c the gap vanishes, and at β c the Hamiltonian is gapless. The way the gap vanishes determines the dynamical critical exponent z of the associated KIM. • The ground state (1) fulfills the area law and does not exhibit any special behavior at criticality. It can be exactly represented as a MPS in 1D, or as a PEPS in higher dimensions [19].
Intrigued by the fact that the states (1) always fulfill the area law, we have attempted to look more closely into their properties, and the properties of their parent HamiltoniansĤ. One possible way is to study the entanglement properties of the excited states ofĤ, but we have chosen another approach. We have generalized the KIMs to quantum models by defining a new class of QMEs, as explained in the previous subsection. These QMEs define new classes of parent Hamiltonians, and the grounds states of these Hamiltonians are expected to be again well described by matrix product or pair entangled projected states. This paper is devoted to the discussion of the critical properties of such Hamiltonians, as well as entanglement properties of their low energy states.

Plan of the paper
The paper is organized as follows. In Secion 2 we remind the readers of basics of KIMs with detailed balance, associated Hamiltonians, and more. In Section 3 we present the main result of this paper: we describe the generalization of KIMs to QMEs, and ‡ We would like to warn our readers against a possible misunderstanding: At issue are irreversible processes whose quantum behavior in general involves mixed states described by density operators. The generators in the pertinent master equations are non-Hermitian operators, precisely due to the irreversible dynamics. However, given detailed balance a master equation can be rewritten such that the generator becomes Hermitian and is then called a "Hamiltonian". Following widespread practice we shall indulge in calling a so transformed master equation a "Schrödinger equation in imaginary time" and talking about eigenstates of the "Hamiltonian"; the resulting formalism then even allows to deal with superposition states, Eq. (1) being the first example. Nevertheless, this paper is about density operators and damped motion. The "state" in Eq. (1) represents the canonical density operator.
we discuss the properties of the QMEs. In particular, we show how the equations for 2 N × 2 N matrix elements of the density matrix split into 2 N equations of the KIMtype for functions of 2 N spin configurations. Detailed balance allows to transform these quasi-KIM equations to a Hamiltonian form, and we derive here the new classes of Hamiltonians. Properties of these Hamiltonians and their low energy states are discussed in Section 4. Here we show that Hamiltonians describing the evolution of the off-diagonal density matrix elements (coherences) are strictly positive §, which implies that the corresponding coherences decay to zero. At criticality we observe effects of critical slowing down of coherences; the Hamiltonians become gapless, implying that at criticality the decay is of the form of an exponential decay times algebraic tails. In some cases, criticality may even lead to survival of coherences for infinite times. In Section 4 we also discuss entanglement properties of the ground and excited states of the corresponding Hamiltonians. Our conclusions and outlook are contained in Section 5.

Basics of kinetic Ising models
As we discuss in Section 4, the methods developed in this paper can be straightforwardly generalized and applied to models other than Ising models: Potts models, classical clock models, or models employing commuting stabilizer operators. For simplicity and concreteness we will limit the discussion to Ising models, with Ising variables σ i = ±1, described quantum mechanically by the commuting Pauli matrices σ z i . We will consider systems with a classical Hamiltonian H(σ) following Markovian dynamics toward the thermal equilibrium state. The dynamics for a Markovian stochastic process is most conveniently formulated for the conditional probability P (σ, t|σ 0 , 0) for the configuration σ at time t, provided the initial configuration was σ 0 at t = 0. This conditional probability allows to calculate all multi-time correlation functions of the process. Note that P (σ, 0|σ 0 , 0) = δ σσ0 , i.e. at the initial time conditional probability is obviously given by the Kronecker delta of the configurations σ and σ 0 . In the following, in order to avoid too many arguments we will consider dynamics of probability of configurations P (σ, t), which fulfills the same equation as P (σ, t|σ 0 , 0), but with a more general initial condition.

Kinetic Ising models with detailed balance
Here, we define in a more detailed way the classical kinetic Ising models and recall some of the literature results concerning particular examples of such models. In general, kinetic Ising models are defined by specifying the so called kinetic (master) equation for probability P (σ, t) of the forṁ where the sum runs over all possible configurations σ. The function w(σ → σ ′ ), hereafter also called transition probability, stands for the probability per unit time that the system changes its configuration from σ to σ ′ . Usually one imposes the so called detailed balance condition (DBC) that says that at equilibrium the probability per unit time of a transition from σ to σ ′ is the same as the probability of transition in the opposite direction, σ ′ → σ. Mathematically, it reads where P eq (σ) is the equilibrium probability distribution P eq (σ) = P (σ, t → ∞).
To support the above general remarks with some more detailed investigations let us now discuss some previously studied examples of such kinetic models. At the beginning let us focus on the one-dimensional Ising spin system, i.e., a chain of N spins uniformly distributed on a line. In this case σ denotes one of the 2 N possible configurations of N spins and can be represented as a N -dimensional vector We restrict our attention to the case in which the behavior of ith spin is local, i.e., depends only on the nearest neighbors (generalizations to local models with next nearest neighbors are straightforward). We also assume that the probability distribution at equilibrium is with H denoting the classical (local) Ising Hamiltonian. In particular, we consider here the ferromagnetic Ising model which in 1D corresponds to In this case the partition function has the explicit form Z N = 2 N (cosh N βJ + sinh N βJ). The simplest possible process that may occur here is a single flip of the ith spin. Schematically this process can be stated as σ → D i σ, where D i denotes the flip at ith position, D i σ = (σ 1 , . . . , −σ i , . . . , σ N ). Also, let w(σ → D i σ) denote the transition probability for that process. The only processes that lead to the configuration σ, appearing on the left-hand side of Eq. (2) are spin flips D i σ → σ for any i = 1, . . . , N . The inverse type of processes can drag the system away from σ. This means that the general kinetic equation (2) can be reduced in this case to a much simpler forṁ The most general form of w(σ → D i σ) in the case that the interaction with both nearest neighbors is symmetric and leads the system to the equilibrium state (4) was shown to be [27]: where γ = tanh 2βJ (notice that the value γ = 0 corresponds to infinite temperature, while γ = 1 to zero temperature), |δ| ≤ 1 and 0 < Γ < ∞. The case of δ = 0 was thoroughly investigated by Glauber, and shown to be solvable in the sense that all the relevant physical quantities can be computed analytically. In particular, the non-equilibrium expectation values and equilibrium correlation functions were determined. Moreover, this model was shown to have the dynamical exponent z = 2. Later on, the above model was treated in a series of This means that the time dependent spin-spin correlation function decay on time scale as t dec behaves as t dec ∝ ξ z , where z is the dynamical critical exponent, and ξ denotes the correlation length; ξ scales as (1 − γ) −1/2 when γ → 1.
KIMs have also been studied in two and higher dimensions, although in these cases there is no known analytical solution. However, using efficient techniques precise numerical results have been obtained for relevant quantities such as the critical dynamical exponent [33].

Associated quantum Hamiltonians
Interestingly, as it was show in e.g. Refs. [27,29,30] the detailed balance condition (3) allows to rewrite the master equation (2) as a Schrödinger equation. In order to see that let us introduce the function φ(σ, t) related to the probability P (σ, t) through P (σ, t) = P eq (σ) φ(σ, t), where as above P eq is the equilibrium distribution. Inserting the latter in the master equation (2) and reorganizing slightly some terms, we arrive atφ with H denoting the real matrix with elements given by Due to the detailed balance condition (3) one sees that H is symmetric and hence Eq. (8) gives the aforementioned Schrödinger equation. Diagonalization of the Hamiltonian H is equivalent to the full solution of the corresponding master equation (2). Many of the previously discussed systems were investigated from this point of view. Below we recall some of the known results for the one-dimensional Glauber model. The Hamiltonian associated to Glauber's master equation (8) with spin rates given by (7) has the form where γ, δ, and Γ are specified as previously; σ z and σ x are the standard Pauli matrices and For δ = 0 it was shown in Ref. [29] (see also Ref. [35] for another approach) that the standard procedure consisting of a Jordan-Wigner transformation followed by Fourier and Bogoliubov-Valatin transformations [37,38] brings the Hamiltonian H(γ, 0) to its diagonal form with non-interacting fermions. The eigenvalues are given by where q k denotes the ordered chain (q 1 < q 2 < . . . < q k ) with each q i from {±π/N, ±3π/N, . . . , ±(N − 1)π/N } for even N and {0, ±2π/N, ±4π/N, . . . , π} for odd N . The ground state, which is here the zero-energy state, is the one given by Eq.
(1). Moreover, the first nonzero eigenvalue is 1 − γ and goes to zero for γ → 1 (zero temperature limit). The fermionic representation allows to go beyond Glauber's seminal work and solve exactly a whole class of master equations that are associated to the Hamiltonians given by Eq. (10). The inverse approach has also been explored: Given a quantum Hamiltonian that can be solved (at least partially), what are the corresponding classical master equations, and do they represent interesting physical systems [39]?

Quantum kinetic Ising models
Now we are in a position to proceed with the quantum generalization of the classical kinetic equations (2). Here we only discuss such a generalization for the Glauber master equation (6) with single-spin flips. However, further generalizations are obviously possible and are left for future work (see also Sec. 5). We start by defining our notation. Naturally, as the computational basis in (C 2 ) ⊗N we take the eigenstates of N -fold tensor product of σ z . The basis consists of 2 N vectors hereafter denoted by |σ ≡ |σ 1 , . . . , σ N (σ = 0, . . . , 2 N − 1) with σ i denoting an eigenvalue of σ z i . After appropriate rescaling we may look at (σ 1 , . . . , σ N ) as the binary representation of the decimal number σ.
Let us consider the following master equation Here Notice that Eq. (13) can be written as with the Lindblad operators given by L i = σ x i [w i (σ z )] 1/2 . It is then clear that we deal only with the dissipative part of the general master equation describing Markov processes.
Nevertheless, further generalizations to the full master equation are possible (see e.g. Refs. [21] and [17]) as, for instance, for the Ising Hamiltonian, the part governed by the Hamiltonian vanishes and such an equation would also reproduce the kinetic equations (2).
In what follows we consider only the spin rates that satisfy detailed balance, i.e., ¶ For the sakes of clarity and simplicity in the case of quantum spin rates we change sightly the notation from w(σ → D i σ) to w i (σ z ). According to this convention w(D i σ → σ) is replaced with w i (D i σ z ).
Let us now proceed with solving the above master equation. As we will see below, our method relies -as in the case of classical kinetic equations -on the observation that the whole equation can be brought to a set of 2 N Schrödinger equations. For this purpose it is convenient to use the isomorphism between matrices from M d (C) and vectors from C d ⊗ C d . More precisely, we can represent the density matrix ̺(t) as This form of ̺(t) allows to rewrite the master equation (13) as the following matrix equation It should be emphasized that operators corresponding to the "tilded" and "nontilded" spins act on the "tilded" and "nontilded" kets in the vectors representation of the density matrix ̺(t) (16) (for instance σ x i σ x i |s | s = σ i |s σ x i | s ). It is evident that the matrix appearing on the right-hand side of Eq. (17) is not Hermitian. In order to bring it to Hermitian form we can use the detailed balance condition (15). This suggests the following transformation (see e.g. Ref. [29]) with H denoting the quantum generalization of the Ising Hamiltonian given by Eq. (5). Application of the above to Eq. (17) leads us to where ]. Now, due to the detailed balance condition, it is clear that σ x i and v i (σ z ) commute, and consequently we have a Schrödinger equation |ψ(t) = −H|ψ(t) with Hermitian H. Let us notice that the above procedure, being just a "quantum" generalization of the procedure described in Sec. 2.2, replaces the problem of solving the master equation for N spins with the problem of solving the corresponding Schrödinger equation for 2N spins.
As we will see below the above equation can be further simplified. Namely, identifying operators that commute with H we can split the equation to a group of 2 N Schrödinger equations. For this aim let us notice that H commutes with σ z i σ z i (i = 1, . . . , N ). We can therefore introduce new variables τ i = σ z i σ z i (i = 1, . . . , N ) which are constants of the motion. Then, one sees that the tilded variables can be expressed by τ = (τ 1 , . . . , τ N ) and σ variables as σ z i = τ i σ z i for any i. In other words we have replaced σ and σ by τ and σ, of which τ is conserved. The Hamiltonian H τ then takes the form where τ σ z denotes τ i σ z i (i = 1, . . . , N ). Let us make a comment on the notation. One sees that there are 2 N different configurations of τ -spins. We label them by the natural numbers 0, . . . , 2 N − 1; since we shall repeatedly have to refer to the two fully homogeneous configurations we simply denote these as τ = 0 (all τ -spins up) and τ = 2 N − 1 (all τ -spins down), without further specifying the association of naturals to configurations.
By definition the τ = 0 configuration represents the equal values of σ and σ spins. Consequently, via Eq. (16), the corresponding Schrödinger equation describes the diagonal elements of the master equation (13) and thus recovers the classical kinetic equations (6). In the remaining cases of τ configurations the σ and σ variables have to differ on at least one position implying that the Schrödinger equation related to any τ = 0 describes the set of 2 N off-diagonal elements of ̺(t).
Thus, by identifying N operators commuting with H we brought the solution of the general master equation (13) to the problem of diagonalization of 2 N Hamiltonians, each of dimension 2 N × 2 N .
Let us now concentrate on a particular choice for the transition probabilities w i (σ z ). In what follows we investigate the quantum version of the rates given by (7) (σ i are replaced with σ z i ), that is with Γ, δ, and γ defined as before. Putting this to Eq. (20) and after some algebra we get where and with A(γ, δ) and B(γ, δ) defined as in Eq. (11) and f (x) = (1/2)(1 + x). One sees that for τ = 0 or τ = 2 N − 1 the function f equals one for i = 1, . . . , N and therefore we obtain the Hamiltonian given in Eq. (10). Let us stress again that the Schrödinger equation corresponding to τ = 0 describes the diagonal elements of the master equation (13) and thus gives the classical kinetic equations (6). All the remaining τ s correspond to the off-diagonal elements of ̺(t).
For the interesting special case δ = 0 all Hamiltonians H τ can be diagonalized using the results of Ref. [41]. This is because for any τ = 0, . . . , 2 N − 1, the Jordan-Wigner transformation [36] brings H τ (γ, δ) to Hamiltonians which are quadratic in the fermion operators This, by virtue of the results of [41], means that diagonalization of any of H τ (γ, 0) reduces to the diagonalization of an N × N matrix. Since, in principle, the latter can be performed numerically efficiently, the Hamiltonians (22) are diagonalizable.
For the case of δ = 0 the Jordan-Wigner transformation gives us Hamiltonians containing terms which are quartic in c i . In this case we can try the transformation proposed by Siggia [35]. For this aim we introduce new bond variables They commute when their indices are different (different bonds) whereas for coinciding indices (same bond) they obey the algebra of Pauli matrices (X 2 i = 1, X i Y i = iZ i etc). The bond variables allow to rewrite the Hamiltonian from Eq. (22) as Generally this is the anisotropic Heisenberg Hamiltonian with some external field. In some particular cases such Hamiltonians are analytically diagonalizable (for a pedagogical review on the Heisenberg model and the Bethe ansatz that is used to solve it see Ref. [42]). It seems that the general case δ = 0 defies exact diagonalization. However, in the following subsection we show that for δ = −1, H τ (γ, δ) can be diagonalized analytically for all values of τ .

Properties of the Hamiltonians associated to the QME
Here we study some of the properties of the Hamiltonians H τ . First, we show that all of them are positive and in a majority of cases even strictly positive. We also discuss the cases in which H τ have zero eigenvalues and thus identify possible stationary states of the master equation (13). Then we study fully soluble case of the socalled energy conserving spin flips, where all Hamiltonians H τ τ = 0, . . . , 2 N − 1 are diagonalizable. Finally, we numerically investigate properties of entropy of ground states of our Hamiltonians.

Zero eigenvalues of H τ -stationary states of the evolution
It is interesting to ask which states survive the evolution, that is, which of the off-diagonal elements of ̺(t) do not decay to zero with t → ∞. This can be done by studying the positivity properties of H τ and in particular their eigenstates corresponding to zero eigenvalues. Generically, as we will see below, for any τ except for τ = 0 and τ = 2 N − 1, the Hamiltonians H τ are strictly positive except for two different (γ, δ)-points, namely, δ = γ = 0 and δ = γ = 1. On the other hand, for τ = 0 or τ = 2 N − 1 it is known [29,28] that the corresponding Hamiltonian has zero eigenvalues (two-fold degenerate) for any value of γ and δ. Let us treat the two cases τ = 0 and τ = 2 N − 1 in a more detailed way. For this purpose let us notice that we may write H τ (γ, δ) as H τ (γ, δ) = i H (i) τ (γ, δ) and study positivity of each H (i) τ (γ, δ). On the other hand, for any τ = 0, 2 N − 1 we may divide all such terms into two groups, the one consisting of H (i) τ (γ, δ) for which τ i−1 = τ i+1 and the one for which τ i−1 = −τ i+1 . More precisely, we can write H τ (γ, δ) as In the first case of equal τ s we have One finds that H (i) τ (γ, δ) has four different eigenvalues, each two-fold degenerate. For −1 ≤ δ ≤ 1 their explicit forms are λ 1 = 0, For the chosen parameter region λ is manifestly positive. On the other hand, to prove nonnegativity of λ In the case of τ i−1 = −τ i+1 we have Here one finds that since vanishes. Without any loss of generality let us assume that f (τ i τ i+1 ) = 0. Then, obviously f (τ i−1 τ i ) = 1 and the above Hamiltonian can be brought to It has two different eigenvalues (each two-fold degenerate) of the form λ ( =) One sees that obviously λ ( =) − is nonnegative it suffices to notice that for 0 ≤ γ ≤ 1 and |δ| ≤ 1 the maximal value of the function of γ and δ appearing under the sign of square root is four and is this value is attained for γ = δ = 0 and δ = γ = 1. This also means that γ = δ = 0 and γ = δ = 1 are the only points for which H (i, =) τ (γ, δ) can have zero eigenvalues. In conclusion, it follows from the above analysis that all H (i) τ (γ, δ)s are positive and therefore the Hamiltonian H τ (γ, δ) is positive for any τ = 0, . . . , 2 N − 1, and parameter region specified by the conditions −1 ≤ δ ≤ 1, 0 ≤ γ ≤ 1, and Γ > 0. Moreover, it follows that for all τ = 0, 2 N − 1, the Hamiltonians H τ (γ, δ) are in general strictly positive except the points δ = γ = 0 and δ = γ = 1. Thus these two points together with two values of τ = 0, 2 N − 1 are the only possible cases when H τ (γ, δ) can have zero eigenvalues. Let us discuss shortly each of these cases.
For τ = 0 and τ = 2 N − 1, as previously noticed, the corresponding Hamiltonian (10) has a zero-eigenvalue eigenstate given by Eq. (1). The case τ = 0 (τ = 2 N − 1) corresponds to the diagonal (anti-diagonal) elements of ̺(t), which as it follows from Eqs. (1) and (18) are of the form (taking into account the normalization) For the zero temperature limit (β → ∞), the above becomes the well-known N -partite Schrödinger cat state (or Greenberger-Horne-Zeilinger state) with |↑ and |↓ denoting the eigenvectors of σ z corresponding to the positive and negative eigenvalue, respectively. Moreover, for γ = 1 the two-fold degeneracy appears in the ground state of H 1 (1, δ) and the second zero-energy eigenstate is Since generically the remaining off-diagonal elements of ̺(t) vanish in the t → ∞ limit (except for the already mentioned values of γ and δ), we have an example of a state that is a stationary state of the master equation and becomes genuine multipartite entangled state in the limit of zero temperature.
The case of δ = γ = 1 is a little bit more difficult. Now, from Eq. (22) it follows that It is clear that this Hamiltonian is diagonal in the standard basis |σ in (C 2 ) ⊗N and thus we can look for the eigenstates among the standard basis in (C 2 ) ⊗N . Interestingly, using the previously introduced bond variables Z i = σ z i−1 σ z i , in the case of periodic boundary conditions this Hamiltonian can be brought to the antiferromagnetic Ising Hamiltonian with magnetic field Let us concentrate on the zero-energy eigenstates of H τ (1, 1). The latter can be found by solving the corresponding equation for eigenvalues of σ z . This, however, due to the fact that in this equation each term under the sum is positive, means solving of the following set of equations (39) where σ i stands for the eigenvalue of σ z i . For instance, for all configurations of τ that consist of blocks of length greater or equal two separated by the domain walls, one of the possible solutions is given by σ i = τ i . This is because in such case the above set becomes f (τ . . . , N ). The only possible triples (τ i−1 , τ i , τ i+1 ) that can appear in the discussed case are ↑↑↑, ↓↓↓, ↑↓↓, ↓↑↑, ↑↑↓, and ↓↓↑. It is clear that for all of them these equations are satisfied. One may also easily check that generically the zero-energy eigenstates are degenerated. As a result it seems that in the case of γ = δ = 1 there is a variety of states ̺(t) that are stationary states of our master equation.

The case of energy conserving spin flips
Here we consider the case of energy-conserving spin flips, that is, when δ = −1. One sees now that the Hamiltonians H τ simplify significantly and read where f τ i = f (τ i−1 τ i+1 ). For τ = 0 and τ = 2 N − 1 this is the isotropic ferromagnetic Heisenberg model with spectrum shifted by ΓN . In the general case, however, H τ (γ, −1) depends on the numbers f i which are either zero or one depending on the configuration τ . To deal with this it suffices to notice that for any configuration of τ = 0, 2 N − 1 some set of indices {i 1 , . . . , i k } exists for which f i k = 0 and between these zeros the f i are constant and equal to one. It is clear from Eq. (40) that such zeros divide the Hamiltonian into a sum of "smaller" commuting Hamiltonians. More precisely if f ij = 0 for j = 1, . . . , k then where S n = [X n , Y n , Z n ]. It follows that S n−1 ·S n commutes with S m−1 ·S m whenever |n − m| ≥ 2. This means that for a given configuration of τ we can split H τ (γ, −1) into a group of commuting isotropic ferromagnetic Heisenberg Hamiltonians with free ends. Such Hamiltonians can be treated using the so-called Bethe ansatz [43]. To visualize what we have just said let us consider the following illustrative example. Let us assume the periodic boundary conditions in Eq. (40) and let the τ configuration together with the corresponding chain f τ i be given by τ ↑ ↑ ↑ ↓ ↓ ↑ ↑ ↑ ↓ ↑ ↓ ↑ ↑ ↑ f i 1 1 0 0 0 0 1 0 1 1 1 0 1 1 (42) Then the Hamiltonian becomes (forgetting about the constant part and putting Γ = 1)

Entropy of the ground state of the off-diagonal Hamiltonians
Here we take a brief detour and study the entanglement between parts of the ground state of the Hamiltonians H τ (γ, δ) that control the dynamics of the off-diagonal terms of the QME. In order to study our system for arbitrary values of δ, we take advantage of the matrix product state structure that the ground state of these Hamiltonians has. For this, we extended the time evolving block decimation (TEBD) algorithm of Vidal [44] to include next-nearest neighbor interaction. We then performed a variable step imaginary time evolution to find the ground state of open boundary chains of length N large enough so that the results become independent of size. Our goal is to compute the bipartite entropy S = tr ρ L log 2 ρ L , where ρ L is the reduced density matrix obtained after tracing out N − L neighboring spins of the chain. For this, instead of obtaining and diagonalizing ρ L , we make use of the Schmidt coefficients that appear explicitly in the TEBD representation of the state + . The system is a chain of 90 spins. Notice how at the point of the flipped τ i 's the bipartite entropy is reduced but does not go to zero. In thin lines, for comparison, is the same entropy but for a configuration with all τ i 's equal to one (the classical KIM, here the diagonal part of the QME).
We study two representative cases of τ configurations, one in which only two τ i 's are different than the rest, and one where half neighboring τ i 's are equal between them and different than the other half -a case that we call a "domain wall". We also concentrate on two relevant values of δ, first δ = 0 (the Glauber model), shown in Figs.
We observe that the maximum value of entropy grows with γ and approaches unity for criticality (γ = 1). The block length L at which the entropy saturates also appears to be rather small, about 5 sites, although we expect this to grow near the critical point.
In all cases of δ a domain wall appears to de-entangle the two parts of the chain, + Any pure state |ψ of a quantum system partitioned into two parts A and B can be written in its Schmidt decomposition form, |ψ =  Figure 3. The same system as Fig. 1, but for δ = γ/(2 − γ), where the classical KIM shows an anomalous dynamical exponent z = 4. In this case, the flipped τ i induces a spike in entropy near the domain wall, indicating some interesting quantum correlation existing between domain walls. Still the curves, in ascending order, are for γ = 0.1, 0.3, 0.5, 0.7, and 0.9. In thin lines, for comparison, is the same system but for a configuration with all τ i 's equal to one (the classical KIM, here the diagonal part of the QME).
giving zero entropy for a partition right at the domain boundary. On the other hand, a single flipped τ i shows some residual entropy. Therefore, the separation of the Hamiltonian into commuting parts shown above for δ = −1 is not possible for general τ configurations.  Figure 4. The same system as Fig. 3, but for a domain wall configuration for τ . In contrast to the case shown in Fig. 2, here the entropy grows sharply before the end of the domain. Nevertheless, it goes to zero exactly at the domain wall, again indicating that in this case the domains are not entangled.
The case with an anomalous dynamical exponent z = 4, δ = γ/(2 − γ), shows some interesting behavior of the entropy near the domain walls or the flipped τ i . In particular, we observe an entropy spike before the end of the domain. The effect is acutely pronounced near criticality.

Conclusions and Outlook
We have presented here a novel class of QMEs that have the following properties (cf. e.g. Refs. [21] and [17]): • The diagonal elements of the density matrix in the "computational" basis follow dynamics of a certain classical kinetic model.
• The dynamics of the off-diagonal matrix elements splits into blocks, described by "kinetic-like" models.
• For models fulfilling DBC, the dynamics can be transformed into a Hamiltonian dynamics in imaginary time.
• The ground and low excited states of the resulting Hamiltonians fulfill area law, and can be well described by MPS or PEPS methods.
Our results suggest several directions of investigations, which we would like to follow in the future: • In the present paper we have focused on Ising spins, and on generalized QME associated with KIMs. Generalizations to models associated with kinetics of more general set of commuting operators, such as stabilizer operators, are possible and interesting. • Several presented models admit exact solutions via Wigner-Jordan transformation or Bethe ansatzá la [35,29], and/or approximate treatment using variational methodsá la [28]. These methods should allow for more rigorous analysis of the novel types of Hamiltonians.
• Especially interesting are the two-spin-flip models such as those that conserve the magnetization [40] or energy [30]. In particular, for the energy conserving model in 1D the generalized QME reads: The high degeneracy of the ground state in this model allows to expect the appearance of long-living coherences. • Generalization to models that do not fulfill DBC such as exclusion models (see Refs. [45,46]) is possible. An example of somewhat "hidden" DBC is the QME of the form where −1 ≤ α ≤ 1 is a free parameter. This model corresponds to 1D anisotropic Heisenberg model.
• Last, but not least, physical implementations of the considered models with ultracold atoms in optical lattices, or ions in trap arrays, or Rydberg atoms are feasible, and will be studied.