Boundary dissipative spin chains with partial solvability inherited from system Hamiltonians
Chihiro Matsui Naoto Tsuji

Optimal Regularity for Fully Nonlinear Nonlocal Equations with Unbounded Source Terms
    Disson S. dos Prazeres and Makson S. Santos
    Received 16 July 2024; accepted 04 September 2024
=====================================================================================





Abstract


Partial solvability plays an important role in the context of statistical mechanics, since it has turned out to be closely related to the emergence of quantum many-body scar states, i.e., exceptional energy eigenstates which do not obey 

the strong version of the eigenstate themalization hypothesis. 
We show that partial solvability of a quantum many-body system can be maintained even when the system is coupled to boundary dissipators under certain conditions. 
We propose two mechanisms that support partially solvable structures in boundary dissipative systems: The first one is based on the restricted spectrum generating algebra, while the second one is based on the Hilbert space fragmentation. 
From these structures,
we derive exact eigenmodes of the Gorini-Kossakowski-Sudarshan-Lindblad equation
for a family of quantum spin chain models with boundary dissipators,
where we find various intriguing phenomena arising from the partial solvability of the open quantum systems, including persistent oscillations (quantum synchronization) and the existence of the matrix product operator symmetry.
We discuss how the presence of solvable eigenmodes affects long-time behaviors of observables in boundary dissipative spin chains based on numerical simulations using the quantum trajectory method.


0
We consider two kinds of partially solvable Hamiltonians coupled to boundary dissipators, for each of which we introduce a mechanism inducing partial solvability under the time evolution given by the Gorini-Kossakowski-Sudarshan-Lindblad equation. 

The first type consists of the Hamiltonians which admit the restricted spectrum generating algebra. We show that several solvable energy eigenstates of such Hamiltonians do not feel the effect of boundary dissipators and therefore become the solvable eigenmodes of the Liouvillians. The corresponding Liouvillians exhibit long-lived oscillations in the expectation values of observables, since they inherit the restricted spectrum generating algebras of the Hamiltonians. 
The second type is provided by the Hamiltonian with the Hilbert space fragmentation with embedded integrability. We show that the boundary dissipators coupled to such a Hamiltonian can be regarded as partial integrability non-violating perturbations in the thermofield double formalism, when they act as quasiparticle baths which dope and absorb quasiparticles. In the solvable subspace, this Liouvillian can be mapped to two decoupled integrable spin chains whose eigenvalues and eigenvectors are exactly calculated via the Bethe ansatz method. 





§ INTRODUCTION

Integrability of isolated quantum systems has been studied for a long time. Since the achievement by H. Bethe <cit.> who has derived the exact eigenfunctions for the Heisenberg spin chain, the method, now called the (coordinate) Bethe ansatz, has become a powerful tool for systematically constructing eigenfunctions of interacting many-body systems.


What lies behind the Bethe ansatz is the decomposability of a many-body scattering into “consistent" two-body scatterings, that is, physics does not depend on a way to decompose a many-body scattering into two-body ones. This property, which is often referred to as the definition of “quantum integrability", is guaranteed by the existence of the R-matrix that solves the Yang-Baxter equation (YBE) <cit.>. 











Once the mathematical background of integrable systems has been understood, various methods have been invented for calculating energy spectra, form factors, correlation functions, etc. The methods include, e.g., the algebraic Bethe ansatz <cit.>, the vertex-operator approach <cit.>, the separation of variables <cit.>, and the off-diagonal Bethe ansatz <cit.>. 

On the other hand, there exists another class of solvable systems, in which only a part of the energy spectrum is analytically accessible. We shall call such a system “partially solvable".

Partial solvability is mostly understood as an extra symmetry of the Hamiltonian that emerges only in a subspace W of the entire Hilbert space ℋ. In other words, the Hamiltonian restricted in W satisfies extra commutation relations induced by the extra symmetry. 

If the restricted Hamiltonian is integrable, 

such symmetry-induced commutation relations are derived from

infinitely many commuting transfer matrices originating from the YBE, leading to the existence of infinitely many conserved quantities. However, in most of the known cases of partially solvable systems, the extra symmetry is specified by a much 
simpler algebraic relation such as “the restricted spectrum generating algebra (rSGA)" <cit.>,



which holds in the subspace W. 
Examples of the partially solvable systems that exhibit the rSGA include the perturbed spin-1 XY model <cit.> and the Affleck-Kennedy-Lieb-Tasaki (AKLT) model <cit.>. 
One can construct energy eigenstates in the solvable subspace
by applying the spectrum generating operator Q (satisfying the rSGA with the Hamiltonian)  to a simply constructed energy eigenstate.



Another known mechanism which may induce partial solvability is the Hilbert space fragmentation (HSF) <cit.>. 
The HSF is defined as the fragmentation of the total Hilbert space into exponentially many invariant subspaces ℋ_r (“Krylov subspaces") induced by the Hamiltonian,



    ℋ = ⊕_r ℋ_r,

provided that such a decomposition is not caused by an obvious local symmetry of the Hamiltonian. 






The fragmentation structure is now mathematically understood in terms of “the commutant algebra”<cit.> for the Hamiltonian, 

which also tells the number of Krylov subspaces induced by the Hamiltonian. 
Although the HSF is not necessary related to the notion of integrability, some models admit both the HSF and integrability (Fig. <ref>(a)) <cit.>. A representative example is the XXC model <cit.>, which has been introduced as a new type of integrable systems without referring to the HSF structure. 
As will be shown in <cit.>, 

the Hamiltonian of the XXC model has the fragmentation of the Hilbert space into exponentially many invariant subspaces, according to “frozen" partial spin configurations (which we call “irreducible strings (IS)" in the main text.)

Surprisingly, the integrable Hamiltonian with the HSF structure can be deformed in such a way that keeps its integrability in a selected subspace among exponentially many of those (Fig. <ref>(b), (c)) <cit.>. The key idea is to choose a perturbation that vanishes

in a selected subspace and hence does not violate 

integrability in this subspace (Table <ref>), although integrablity in the entire Hilbert space is in general lost by such a perturbation.

Based on a similar idea, we will introduce a variety of perturbations that keeps integrability in a selected subspace in the main text. A remarkable fact is that partial solvability in a selected subspace holds even under site-dependent integrability-breaking perturbations, since these perturbations are irrelevant in a selected subspace (a similar idea for Hamiltonians with rSGA can be found in <cit.>). 

0

0




Partially solvable systems are now intensively studied especially in the context of thermalization. 

For instance, the solvable subspace W is an invariant subspace of the Hamiltonian, and therefore, any state in W never reaches the other subspaces during time evolution. This is a typical example of “weak ergodicity breaking" in the Hilbert space, which may be considered as a necessary condition for the emergence of “quantum many-body scar (QMBS) states" <cit.>, i.e., non-thermal states in a non-integrable 

system. Indeed, many of the QMBS states have been found to be exactly solvable energy eigenstates of non-integrable systems. Several examples can be found in <cit.>. 
Another remarkable feature in partially solvable systems is a persistent oscillation of local observables <cit.>. This phenomenon is understood as a consequence of a large overlap between an initial state and solvable energy eigenstates forming equally-spaced energy spectra imposed by the rSGA.

Existence of long-lived oscillations implies that the system neither thermalizes nor relaxes

to any steady state. 

Motivated by these atypical behaviors of partially solvable systems, we focus on a question of whether an open quantum system can also admit partial solvability, and if so, what are characteristic phenomena in partially solvable open quantum systems. Let us consider open quantum systems that evolve according to the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation <cit.>. 






The GKSL equation constitutes the most general form of the completely-positive-trace-preserving map under the assumptions that time evolution is Markovian (sometimes one also assumes an initial state to be a product state of the system and environment). 
Recent progress unveils several solvability mechanisms for the GKSL equation, which are classified into two types (Table <ref>): The first one consists of completely solvable Liouvillians whose spectra are fully accessible by analytic methods <cit.>, while the second class consists of partially solvable Liouvillians whose spectra are only partially accessible <cit.>. 
Most of the known partially solvable Liouvillians are solvable only for the steady state, leaving the other eigenmodes including the slowest decaying mode unsolvable (see, e.g., Ref. <cit.> and references therein). Only a few examples are known as partially solvable Liouvillians with a no-less-than-two-dimensional solvable subspace induced by rSGA <cit.>, but in those examples the dissipators are coupled to all the sites of the system.



In this paper, we are especially interested in extending the notion of partial solvability for closed quantum systems to open quantum systems. Our target is an open quantum system with partial solvability inherited from a partially solvable system Hamiltonian, whose solvable energy eigenstates are robust against boundary dissipators.
One way to realize such an open quantum system is to employ a partially solvable system with rSGA, some of whose solvable eigenstates vanish by all the boundary dissipators.




We find that those solvable energy eigenstates may exist for the AKLT-type Hamiltonians when they couple to boundary dissipators injecting spin-2 quasiparticles. 
Another way is to employ a system Hamiltonian with the HSF structure and absorb the boundary dissipators as integrability-preserving

perturbations in a selected subspace. 



We find that there exist site-dependent perturbations that keep integrability in the selected subspace, and hence certain boundary dissipators can be interpreted as integrability-preserving

perturbations in the selected subspace of the integrable doubled spin chain (Fig. <ref>(c)). 

 As we will obtain in the main text, such a system may be realized for the perturbed XXC model whose edges are attached to quasiparticle baths. 
In these partially solvable open quantum systems, the solvable subspace of the system Hamiltonian is 

inherited to a solvable subspace of the Liouvillian. Therefore, the partially solvable Liouvillians constructed in the above ways have not only the solvable steady states but also the solvable eigenmodes. As a result, we observe characteristic behaviors such as persistent oscillations, which can never be obtained for non-integrable Liouvillians. 

The rest of this paper is organized as follows. In Sec. <ref>, we give a review of partial solvability for closed quantum systems, and present our original results on site-dependent perturbations and an integrable subspace encoded by the period-three IS (see Sec. <ref>).

We explain two mechanisms of the partial solvability, the rSGA and HSF, together with some observations about the mechanism to incorporate the HSF and the integrability-preserving

perturbations in a selected subspace. Section <ref> is the main part of this paper, which is devoted to partial solvability of open quantum systems. We show how the partial solvability of the system Hamiltonian can be inherited to open quantum systems, provided that the system shows either the rSGA or HSF. 
In the latter case, we take “the thermofield double (TFD) formalism" <cit.>, which maps a density matrix defined in ℋ⊗ℋ^† to a vector in the doubled Hilbert space ℋ⊗ℋ^*.
The action of the Liouvillian can be expressed as the Hamiltonian for the two decoupled integrable spin chains in the solvable subspace, if the quantum jump terms are irrelevant due to the HSF. 








We also demonstrate numerical results for some of the models to see how solvable eigenmodes in open quantum systems affect long-time behaviors of observables.
The concluding remarks are given in Sec. <ref>, in which some open questions and possible future works are listed.




0
In this paper, we are especially interested in the partially solvable open quantum systems whose partial solvability comes from bulk partial solvability. Our focuses is the partially solvable spin-1 model equipped with the rSGA. Then we set the quasiparticle baths on its ends, in such a way that the system feels quasiparticles are coming in from both ends. 
Under this situation, one may naively expect that the system reaches the fully-polarized steady state and this is the only steady state. But what we have found is that any diagonal ensemble consisting of solvable energy eigenstates can be the steady state. 
Emergence of this unexpected infinitely degenerate steady states is induced by the existence of the exactly solvable subspace whose elements vanish under the dissipators. 
Such density matrices which do not feel the dissipators has first been introduced for the optimal system by the name the dark states <cit.>, since those do not emit photons. Nowadays, the dark state |D ⟩ is defined as a state which satisfies

    H |D ⟩ = E |D ⟩,     
        L_α |D ⟩ = 0,  ∀α,

which are the sufficient conditions for the state |D ⟩ to vanish by the dissipators. We instead propose another mechanism to produce the state which does not feel the dissipators, by showing that the constructed steady states vanish by the dissipators, although they survive by each quantum jump operator. 

It is natural to expect that, by the fact that there exists the solvable subspace which does not feel the dissipators, the elements in this solvable subspace inherits many interesting bulk properties observed as a result of the rSGA structure in the system Hamiltonian. One representative feature is the persistent oscillations, corresponding to the persistent oscillations obtained for the closed quantum systems equipped with the rSGA-based QMBS <cit.>. Actually, similar phenomena have been first reported for another partially solvable quantum system <cit.>, although its partial solvability does not come from the bulk partial solvability but from the cooperation between the Hamiltonian and dissipation operator. 
Consequently, we have theoretically confirm that the persistent oscillations emerge for our open system only when the initial state is prepared as a superposition of solvable energy eigenstates of H, contrary to the known example <cit.> in which the similar oscillation can be observed by stating from any superposition of energy eigenstates with real coefficients. Besides, we numerically verify the long-lived oscillations for the ℤ_2-invariant initial state, which indeed is the inherited feature from the bulk <cit.>. 

Besides the non-integrable models, several integrable models can also admit the HSF structure <cit.>. Curiously, some models succeeded to keep their integrability even in several Krylov subspaces even by adding integrability-breaking perturbations <cit.>. Suppose we look for the perturbations that keeps integrability in a given subspace ℋ_0. Then these perturbations can be any as long as they vanish under the projector onto the integrable subspace ℋ_0. That means, these integrability-breaking perturbations can be site-dependent, unlike the conventionally considered examples most of which are translationally invariant, even by admitting several choices for impurities. 
Therefore, our plan for this paper is first to show the construction of integrable models with the HSF, and then provide the examples of partially integrability-breaking perturbations for each model. 

Another interesting aspect of the partially solvable models with site-dependent perturbations can be observed in the relation to open quantum systems. Especially, recent developments in the open quantum systems whose time evolution obeys the Gorini2013Kossakowski2013Sudarshan2013Lindblad (GKSL) equation achieved classification of their solvability <cit.>. 
In spite of great success in clarifying their solvability mechanism, translation invariance is often imposed when discussing solvability of Liouvillians, i.e. most of the considered examples couple to dissipators at every site. Such a setup is not very useful from the viewpoint of experimental accessibility. 
For this reason, we aim to realize the partially solvable open quantum system coupled to dissipators only at its boundaries, which we expect is easily examined experimentally. 

In order to examine time evolution of open quantum systems which obeys the GKSL equation, “the thermofield double (TFD) vector" formalism <cit.> is often employed, in which a density matrix defined in the ℋ⊗ℋ^† is mapped to a vector in the doubled Hilbert space ℋ⊗ℋ^* by an isomorphism 

    φ:    ℋ⊗ℋ^†→ℋ⊗ℋ^*, 
       ρ = ∑_m,n p_m,n |m ⟩⟨ n| 
        ↦ |ρ⟩⟩ = ∑_m,n p_m,n |m ⟩⊗ |n ⟩^*.

Subsequently, the Liouvillian is mapped to the non-Hermitian Hamiltonian in the doubled Hilbert space. 
In the TFD formalism, a Liouviillian with boundary dissipators is mapped to a periodic non-Hermitian Hamiltonian with two impurities corresponding to boundary dissipators. Thus, one can regard this boundary-driven open system as the closed periodic chain with local impurities. 
This fact evokes us to apply the idea of site-independent partially integrability-breaking perturbations arises in the context of the HSF to constructing partially solvable Liouvillians. 
We show that, indeed, several choices of boundary dissipators can match the form of partially integrability-breaking perturbations. 

In this paper, we propose two mechanisms for robustness of the invariant subspace consisting of QMBS of the bulk Hamiltonian against certain boundary dissipators. The first one is “the boundary cancellation mechanism" on QMBS in the matrix product forms. The second one is the HSF with the invariant subspace in which the dissipators are irrelevant. 

Quantum thermalization is one of the most mysterious phenomena in statistical mechanics. People have devoted efforts to understand quantum thermalization based on quantum mechanics over decades and recently big progresses have been achieved by the reinterpretation of eigenstate thermalization hypothesis (ETH) <cit.>. 

While classical thermalization is understood based on two hypotherical mechanisms the principle of equal probability, which says every microscopic states in the energy shell are realized with the equal probability, and typicality, which requires for almost all microscopic states in the energy shell to be macroscopically undistinguishable from the thermal state, quantum thermalization is now believed to be understood by the following two mechanisms. The first one is ergodicity in the Hilbert space, which corresponds to the principle of equal probability for classical statistical mechanics. This guarantees for a quantum system to relax to a certain steady state. As fluctuations of a local observable by time around its time average are evaluated by 

    lim_T →∞1/T∫_0^T |⟨ O ⟩(t) - ⟨ O ⟩_ d |^2  dt ≤|O||^2/D_ eff,     
        D_ eff = ( ∑_j |c_j|^4 )^-1

the right-hand side, the upper bound of fluctuations, indeed vanishes in the thermodynamic limit if the initial state is prepared to be a superposition of many enough energy eigenstates <cit.>, which guarantees ergodicity in the Hilbert space. 
The second one is the eigenstate thermalization hypothesis (ETH) <cit.>, which is a quantum analog of classical typicality. ETH states that all or almost all, depending on how much strong the statement is, energy eigenstates cannot be distinguished from the thermal state, as long as their eigenenergies are in a macroscopically small interval. This is rather strong statement, and indeed, recent studies show that there exist exceptional non-thermal energy eigenstates in various systems <cit.>. 

These exceptional non-thermal energy eigenstates usually appear in systems which do not thermalize, such as integrable sytems or many-body localized systems <cit.>. However, recently it has been discovered that such non-thermal energy eigenstates show up in systems which do thermalize <cit.>. These states are called quantum many-body scars (QMBS), which are intensively studied these five years. 

The first paper which reported the emergence of QMBS is <cit.>, in which persistent revivals have been observed in the domain-wall density on the Rydberg atom quantum chain. Nowadays, emergence of such long-lived revivals is theoretically understood from the algebraic viewpoint as a result of the hidden restricted spectrum generating algebra (rSGA) <cit.>. The restricted spectrum generating algebra is a kind of dynamical symmetry which requires for the presence of a local operator Q which does not commute with the Hamiltonian but satisfies

    [H,  Q] - ω Q |_W = 0,   ω∈ℝ

in a subspace of the entire Hilbert space W ⊂ℋ. 
Therefore, once we find an exactly solvable energy eigenstate ψ_0 ⟩, we can subsequently construct the other energy eigenstates from ψ_0 ⟩ by applying Q. In other words, the solvable energy eigenstates constructed in this way form the invariant subspace W. These energy eigenstates exhibit the spectrum equally spaced by ω, which perfectly explains why the oscillations we just mentioned emerge. 
Since the solvable subspace constructed in this way is negligibly small in the thermodynamic limit, the model equipped with the rSGA weakly violates ergodicity in the Hilbert space. 
For this reason, exactly solvable energy eigenstates in a partially solvable model are the candidates of QMBS. 

Another characteristic feature of QMBS is relatively small entanglement entropies compared to those of the thermal state, which grows by the volume law. While this has been numerically tested for most cases, some partially solvable models are known to have the exact zero-energy eigenstate in a matrix product form and quasiparticle excitation states on the top of it. As the entanglement entropy of a matrix product state is bounded from above by its bond dimension, these states have constant entanglement entropies as long as the bond dimensions are fixed. 
Indeed, a variety of theoretical models has been proposed as the examples of partially solvable models due to the rSGA <cit.>, all of which have the matrix product zero-energy eigenstate. 

Now we turn to open quantum systems. Stimulated by the recent developments in understanding QMBS, it is natural to ask what are the counterparts to QBMS in open systems.
The time evolution of the open quantum system is governed by the Liouvillian ℒ: 

    d/dtρ(t) = ℒ(ρ(t)).

Actually, solvability of open quantum systems has been intensively studied these days as well, especially when the system obeys the Gorini2013Kossakowski2013Sudarshan2013Lindblad (GKSL) equation. 



Solvability of the GKSL equation is classified into two types.  Liouvillians of the first type are completely solvable by ***. In the second type, on the other hand, systems are partially solvable in which only steady states are exactly solvable. 
In this paper, we look for more general examples of partially solvable open systems in which, besides the steady states, several eigenmodes are exactly solvable. 

A few examples of partially solvable Liouvillians can be found in the literatures <cit.>. The first one discusses ***

The Hilbert space fragmentation (HSF) is known as one of the mechanisms which induces the non-thermalizing property of non-integrable closed quantum systems. 
The HSF has first been introduced in <cit.> and the mechanism for their emergence has been formulated in the context of the models with dipole-momentum conservation <cit.>. 
One can say that the Hamiltonian exhibits the HSF when the Hilbert space are divided into the exponentially many invariant subspaces (“the Krylov subspaces") of the Hamiltonian 

    ℋ = ⊕_r ℋ_r

and this decomposition is not caused by the obvious local symmetry of the Hamiltonian. 
Later in <cit.>, the mathematical structure behind the HSF has been explained by introducing “the commutant algebra", which also tells the number of Krylov subspaces induced by the Hamiltonian. 

Besides the non-integrable models, several integrable models can also admit the HSF structure <cit.>. Curiously, some models succeeded to keep their integrability even in several Krylov subspaces even by adding integrability-breaking perturbations <cit.>. Suppose we look for the perturbations that keeps integrability in a given subspace ℋ_0. Then these perturbations can be any as long as they vanish under the projector onto the integrable subspace ℋ_0. That means, these integrability-breaking perturbations can be site-dependent, unlike the conventionally considered examples most of which are translationally invariant, even by admitting several choices for impurities. 
Therefore, our plan for this paper is first to show the construction of integrable models with the HSF, and then provide the examples of partially integrability-breaking perturbations for each model. 

Another interesting aspect of the partially solvable models with site-dependent perturbations can be observed in the relation to open quantum systems. Especially, recent developments in the open quantum systems whose time evolution obeys the Gorini2013Kossakowski2013Sudarshan2013Lindblad (GKSL) equation achieved classification of their solvability <cit.>. 
In spite of great success in clarifying their solvability mechanism, translation invariance is often imposed when discussing solvability of Liouvillians, i.e. most of the considered examples couple to dissipators at every site. Such a setup is not very useful from the viewpoint of experimental accessibility. 
For this reason, we aim to realize the partially solvable open quantum system coupled to dissipators only at its boundaries, which we expect is easily examined experimentally. 

In order to examine time evolution of open quantum systems which obeys the GKSL equation, “the thermofield double (TFD) vector" formalism <cit.> is often employed, in which a density matrix defined in the ℋ⊗ℋ^† is mapped to a vector in the doubled Hilbert space ℋ⊗ℋ^* by an isomorphism 

    φ: ℋ⊗ℋ^†→ℋ⊗ℋ^*, 
       ρ = ∑_m,n p_m,n |m ⟩⟨ n| 
        ↦ |ρ⟩⟩ = ∑_m,n p_m,n |m ⟩⊗ |n ⟩^*.

Subsequently, the Liouvillian is mapped to the non-Hermitian Hamiltonian in the doubled Hilbert space. 
In the TFD formalism, a Liouviillian with boundary dissipators is mapped to a periodic non-Hermitian Hamiltonian with two impurities corresponding to boundary dissipators. Thus, one can regard this boundary-driven open system as the closed periodic chain with local impurities. 
This fact evokes us to apply the idea of site-independent partially integrability-breaking perturbations arises in the context of the HSF to constructing partially solvable Liouvillians. 
We show that, indeed, several choices of boundary dissipators can match the form of partially integrability-breaking perturbations. 

In this paper, we propose two mechanisms for robustness of the invariant subspace consisting of QMBS of the bulk Hamiltonian against certain boundary dissipators. The first one is “the boundary cancellation mechanism" on QMBS in the matrix product forms. The second one is the HSF with the invariant subspace in which the dissipators are irrelevant. 






§ PARTIALLY SOLVABLE CLOSED SPIN CHAINS



In this section, we mostly give an overview of partial solvability for closed quantum systems. 

In Sec. <ref>, we propose the existence of site-dependent perturbations to a certain integrable Hamiltonian that preserve integrability in a selected subspace, which is our original result. 


We mainly consider

s=1 spin chains with nearest-neighbor interactions as an example. The Hamiltonian can be written as

    H = ∑_j=1^N h_j,j+1,     
        h_j,j+1 = 1⊗⋯⊗j,j+1h⊗⋯⊗1, 
       h = ∑_s,t,s',t'=0^2 h^s,s'_t,t' |tt' ⟩⟨ ss'|,

where h is a local Hamiltonian acting on two neighboring s=1 spins whose states are labeled by |tt'⟩ (t,t'=0,1,2), and N is the number of lattice sites.

Among several mechanisms that produce partial solvability of closed quantum systems, we focus on the restricted spectrum generating algebra (rSGA) and the Hilbert space fragmentation (HSF). 




 §.§ Restricted spectrum generating algebra
 

The notion of the spectrum generating algebra (SGA), or also called the dynamical symmetry, has been introduced in various contexts <cit.>. In this paper, we shall say that the model has the SGA if there exists a spectrum generating operator Q^† that satisfies an algebraic relation,  

    [H, Q^†] - ℰ Q^† = 0,

for a real constant ℰ. 
If a model has the SGA and some of its energy eigenstates are known, one can construct towers of eigenstates by applying the operator Q^† to those known states repeatedly. 

The observed energy spectrum is then equally spaced with the interval ℰ due to the algebraic relation (<ref>). 

One of the simplest examples that show the SGA is a free-fermion model. The Hamiltonian (denoted by H_ FF) is diagonalized in the momentum space,

    H_ FF = ∑_k Λ_k η_k^†η_k,

with real eigenvalues Λ_k and a fermion creation operator η_k^†,
and thus satisfies the SGA,

    [H_ FF, η_k^†] = Λ_k η_k^†,

due to the anti-commutation relations for η_k^† and η_k. 
Then all the energy eigenstates are created by applying the fermion creation operators η_k^† to the obvious vacuum state |0 ⟩, i.e., the zero-energy eigenstate,  

    H_ FF η_k_1^†⋯η_k_n^† |0 ⟩ = (Λ_k_1 + ⋯ + Λ_k_n) η_k_1^†⋯η_k_n^† |0 ⟩.

In this example, the fermion creation operator plays a role of the spectrum generating operator for each mode k. 

Sometimes there appears the SGA only in a subspace W of the entire Hilbert space ℋ,

    [H, Q^†] - ℰ Q^†|_W = 0,    W ⊂ℋ.

We call this type of the SGA structure that emerges only in the subspace W “the restricted spectrum generating algebra (rSGA)". 
As in the case of the SGA that holds in the entire Hilbert space, one can construct a tower of energy eigenstates in the subspace W by repeatedly applying the operator Q^† to an obvious energy eigenstate |Ψ_0 ⟩ (if it may exist).

One of the simplest examples equipped with the rSGA is the perturbed spin-1 XY model <cit.>, whose local Hamiltonian is given by 







    h_j,j+1^XY = J/2 (S_j^+  S_j+1^- + S_j^-  S_j+1^+) + m/2 (S_j^z  + S_j+1^z),

where S_j^±=S_j^x± iS_j^y and S_j^z are spin-1 operators, and J and m are real coefficients.
This model is considered to be non-integrable, and therefore, the energy eigenstates are in general not exactly solvable. However, one can easily find that the fully-polarized state |22… 2 ⟩ (and |00… 0 ⟩) is obviously an energy eigenstate with the eigenenergy -mN (mN). Accordingly, one can derive some of the excited states exactly

in the form of (Q^†)^n |22… 2 ⟩ (that allows a quasiparticle picture),
where Q^† creates quasiparticles of spin-2 magnons carrying momentum k=π


    Q^† = ∑_x=1^N (-1)^x (S_x^+)^2.

The spin-2 magnon creation operator Q^† then satisfies the rSGA,

    [H_XY, Q^†] + 2m Q^†|_W = 0,

in the subspace W

spanned by the quasiparticle excitation states

    W =  span{ (Q^†)^n |22… 2 ⟩}_n ∈{0,1,…,N}.

The rSGA produces the equally-spaced eigenvalue spectrum in the solvable subspace W, which can

explain the persistent oscillation observed in the Loschmidt echo <cit.>. This implies that the perturbed spin-1 XY model never relaxes to any steady state, if the initial state has large enough overlap with the solvable subspace W. 


A more involved example exhibiting the rSGA is a family of the spin-1 Affleck-Kennedy-Lieb-Tasaki (AKLT)-type model <cit.>. The local Hamiltonian is given by 







    h_j,j+1^ AKLT =   1/2 h^00_00 (S_j^x S_j+1^x + S_j^y S_j+1^y + S_j^z  S_j+1^z) 
       - ( 1/2 h^00_00 + a_1^2/a_0 a_2 h^11_11) (S_j^x S_j+1^x + S_j^y S_j+1^y + S_j^z  S_j+1^z)^2 
       - (  1/2 h^00_00 - a_1^2/a_0 a_2(a_1^2/a_0 a_2 - 1 ) h^11_11) (S_j^x S_j+1^x + S_j^y S_j+1^y)^2 
       - ( h^00_00 + ( a_1^2/a_0 a_2 - 1 ) h^11_11) (S_j^z  S_j+1^z)^2 
       + ( 1/2 h^00_00 + ( a_1^4/a_0^2 a_2^2 - 1 ) h^11_11) ((S_j^z)^2 + (S_j+1^z)^2) 
    		+ ( 1 - 2 a_1^4/a_0^2 a_2^2) h^11_11,

in which h^00_00 and h^11_11 are free real parameters, while a_0, a_1, and a_2 are free complex parameters under the condition that a_1^2/(a_0 a_2) ∈ℝ. Note that the original AKLT model is realized by choosing h^11_11 / h^00_00 = 2/3 and a_0 = -√(2) a_1 = -a_2 = √(2/3). 
Although this Hamiltonian is non-integrable, it has been known for a long time that the ground state and some of the excited states are exactly solvable <cit.>, to which a new list of solvable excited states has been added recently in the context of the QMBS states <cit.>. 


The zero-energy state
of the AKLT-type model is written in the form of the matrix product state,

    |Ψ_0 ⟩ = ∑_m_1,… m_N ∈{0,1,2} tr_a (A_m_1⋯ A_m_N) |m_1 … m_N ⟩,

with the frustration-free condition

    h_j,j+1^ AKLTA⃗_j A⃗_j+1 = 0,   
        A⃗_j = [ A_0; A_1; A_2 ]_j,

which holds for j = 1,2,…,N. 

The matrix-valued elements A_0, A_1, and A_2 are given by the Pauli matrices 


    A_0 = a_0 σ^+,    
        A_1 = a_1 σ^Z,   
        A_2 = a_2 σ^-, 
       σ^± = σ^x ± i σ^y,

in which 
the coefficients a_0, a_1, and a_2 are the same as in the Hamiltonian (<ref>). Therefore, the bond dimension of this matrix product state is two. 


Besides the ground state, it has been found that several excited states are solvable, most of which admit the quasiparticle description <cit.>. Here we focus on the excited states expressed by the spin-2 magnons carrying momentum k = π. They are created by the operator given in (<ref>), 



which also satisfies the rSGA for the AKLT-type Hamiltonian,



    [H_ AKLT, Q^†] - 2h^00_00 Q^†|_W = 0,


in the subspace W ⊂ℋ spanned by the quasiparticle excited states 


    W =  span{ (Q^†)^n |Ψ_0 ⟩}_n ∈{0,1,…,⌊ N/2 ⌋}.

Thus, the quasiparticle creation operator plays a role of the spectrum generating operator that creates a tower of solvable energy eigenstates on top of the matrix product 
zero-energy 
state (<ref>). 

Recent studies on partial solvability have provided more formal understanding about the emergence of the rSGA in terms of quasisymmetry <cit.> for the former example and its deformation for the latter example <cit.>. 






 §.§ Hilbert space fragmentation

Hilbert space fragmentation (HSF) is characterized by exponentially many block-diagonal structures of the Hamiltonian H, which are caused by non-obvious symmetries of H. Due to the block-diagonal structure of the Hilbert space, each subsapce is never accessed from the other subspaces by time evolution. 
Among several mechanisms for HSF <cit.>, we focus on the fragmentation observed for the Hilbert space of spin chains due to the presence of “frozen" spin configurations under the action of a Hamiltonian. Those “frozen" spin configurations are called  “irreducible strings" (IS) <cit.>, which are associated with a non-obvious symmetry of the Hamiltonian.

In order to explain the HSF induced by IS, let us 
consider spin chains with arbitrary spin-s
(instead of spin-1).
Suppose that we have a spin-s chain with nearest-neighbor interactions, whose local Hamiltonian is given by

    h^ HSF   = ∑_s,s' ∈ A h^s,s'_s,s' |ss' ⟩⟨ ss'| 
        + ∑_t,t' ∈ B h^t,t'_t,t' |tt' ⟩⟨ tt'| 
          + ∑_s ∈ A, t ∈ B(h^s,t_t,s |ts ⟩⟨ st| + h^t,s_s,t |st ⟩⟨ ts| + h^s,t_s,t |st ⟩⟨ st| + h^t,s_t,s |ts ⟩⟨ ts|).

Here A and B represent subsets of the labels for local states A,B ⊂{0,1,… 2s}, which satisfy A ∪ B = {0,1,… 2s} and A ∩ B = ∅. We also call the labels for local states “species". 
Since the Hamiltonian given in the form of Eq. (<ref>) never exchanges the species in each subset A or B, the configuration (i.e., IS) in each subset A and B is “frozen" under the action of the Hamiltonian.
The existence of such a frozen partial configuration causes fragmentation of the Hilbert space, i.e., the the Hilbert space is fragmented according to the partial configurations in the subsets A and B. 

In each of the fragmented subspaces, one can observe that a spin-1/2 model (with two local states) is embedded in the following way.

The entire Hilbert space of the spin-s chain is ℂ^(2s+1)N, which consists of the tensor product of N local linear spaces ℂ^2s+1 spanned by the (2s+1) basis vectors |0 ⟩, |1 ⟩,…,|2s+1 ⟩ corresponding to the spin degrees of freedom. In this basis, a state in ℂ^(2s+1)N is labeled by the spin configuration. 
Alternatively, one can use another basis 
with local states labeled by A and B, and configurations realized within A and B

(Fig. <ref>). 


Suppose that the subset A consists of N_A species and B consists of N_B (= (2s+1) - N_A) species. The total degree of freedom is then calculated as 

    ∑_n=0^N [ N; n ]· N_A^n N_B^N-n = (N_A +N_B)^N,

which matches the dimension of the Hilbert space in the first description (dim ℂ^(2s+1)N), indicating that the two different bases describe the same Hilbert space. 
For the Hamiltonian with HSF (<ref>), 
the configurations within A and B (i.e., IS) are frozen, while the labels A and B are unconstrained. The (A, B) degrees of freedom form the N-fold tensor product of the two-dimensional local linear space ℂ^2.



The projection onto such a 2^N-dimensional subspace 
is defined by restricting to the configurations in A and B specified by IS (denoted by a projector P_ IS) and then by identifying all the local states in each of the subspaces A and B (Fig. <ref>).





In this way, the Hamiltonian (<ref>) is reduced to a spin-1/2 chain with the nearest-neighbor interactions. 











  §.§.§ Completely integrable case

Although the HSF is not necessarily associated with the notion of solvability, it is often possible to embed solvability in one or a small number of subspaces. 
From now on, let us for simplicity come back to spin-1 systems.

One representative example is the spin-1 XXC model <cit.> with the local Hamiltonian,

    h^XXC = coshη(∑_s,s' ∈{0,2} |ss' ⟩⟨ ss'| + |11 ⟩⟨ 11| )
        + ∑_s ∈{0,2}( |s1 ⟩⟨ 1s| + |1s ⟩⟨ s1| ).

The Hamiltonian exhibits the HSF, without changing the configuration of the labels 0 and 2, and thus the IS for the XXC model is the partial configuration formed by 0 and 2. This local Hamiltonian indeed belongs to the class of (<ref>), by choosing the subsets A = {0,2} and B = {1}. 

Although the projector onto the subspaces specified by the IS generally has a site-dependent complicated expression, sometimes it can be expressed in a simple form.
One example is the projector onto the direct sum of the subspaces encoded by the fully-polarized ISs, … 0000 … and … 2222 …, of length n (n=0,…,N). The projectors onto the corresponding subspace are expressed by simple tensor products of the physical spaces,

    P^(0)_ pol = ⊗_j=1^N (|0 ⟩⟨ 0| + |1 ⟩⟨ 1|)_j,     
        P^(2)_ pol = ⊗_j=1^N (|2 ⟩⟨ 2| + |1 ⟩⟨ 1|)_j,

respectively.
Another example is the projector onto the subspace encoded by the alternating IS, … 0202 … <cit.>,









    P_ alt 
        = ∑_0pt(α_1,β_1),…,(α_N,β_N)∈{0,1,2}^2 tr_a (M_α_1,β_1⋯ M_α_N,β_N)
        |α_1 …α_N ⟩⟨β_1 …β_N|
        - ⊗_j=1^N (|1 ⟩⟨ 1|)_j, 
       M_0,0 = σ_a^+,    M_1,1 = 1_a,    M_2,2 = σ_a^-,    M_α,β = 0    (α≠β),


in which σ_a^± and 1_a represent the Pauli matrices and the two-by-two unit matrix acting non-trivially on the auxiliary space, respectively. Here the trace tr_a is taken only over the auxiliary space.


This projector cannot be written in a simple tensor product form, but instead can be written in an (almost) matrix product form by introducing the two-dimensional auxiliary space. 
Note that the projectors P_ pol^(0), P_ pol^(2), and P_ alt are not orthogonal to each other, but have a common one-dimensional subspace, span{⊗_j=1^N |1 ⟩_j}.

Besides the HSF structure, the XXC model also exhibits integrability, which is guaranteed by the existence of the R-matrix solving the Yang-Baxter equation (YBE), 

    R_1,2(λ_1,λ_2) R_1,3(λ_1,λ_3) R_2,3(λ_2,λ_3) = R_2,3(λ_2,λ_3) R_1,3(λ_1,λ_3) R_1,2(λ_1,λ_2)
           (λ_1,λ_2,λ_3 ∈ℂ),


defined in the three-fold tensor product of the linear spaces V_1 ⊗ V_2 ⊗ V_3. We denote the R-matrix that acts non-trivially on the ith and jth sites by R_i,j, e.g.,

    R_1,2(λ_1,λ_2) = R(λ_1,λ_2) ⊗1_3.


The explicit form of the R-matrix for the XXC model can be found in <cit.>. 

Integrability of the XXC model is inherited by the projected Hamiltonian onto the subspace specified by a certain IS. Regardless of the choice of IS, any projected Hamiltonian onto the subsapce specified by a certain IS results in the same spin-1/2 integrable XXZ Hamiltonian by identifying |0⟩ and |2⟩, 

    h_j,j+1^XXCP_ ISℋ∖{ |0⟩,|2⟩}^N⟼coshη( |↑↑⟩⟨↑↑ | + |↓↓⟩⟨↓↓| )
        + ( |↑↓⟩⟨↓↑| + |↓↑⟩⟨↑↓| ),

where we assign up spins to the sites belonging to the subset A and down spins to the sites belonging to the subset B. 
This also indicates that the spin-1 XXC Hamiltonian can be diagonalized sector by sector labeled by an IS via the spin-1/2 XXZ Hamiltonian, instead of being diagonalized directly via the XXC R-matrix. 





0
Since in the subspace with a fixed IS the XXC Hamiltonian only distinguishes whether each local state belongs to A or B without distinguishing configurations realized in each of the subset A and B, we can consider the projection which maps the XXC Hamiltonian to the reduced Hamiltonian in a given subspace by identifying all configurations in each subsets A or B 

    h_XXCP_ ISℋ∖ (A^⊗ N B^⊗ N)⟼coshη( |AA ⟩⟨ AA| + |BB ⟩⟨ BB| )
        + ( |AB ⟩⟨ BA| + |BA ⟩⟨ AB| ).

This is nothing but the XXZ Hamiltonian with anisotropy coshη. 









Thus, we can diagonalize the spin-s XXC Hamiltonian sector by sector labeled by IS as the spin-1/2 XXZ Hamiltonian, instead by diagonalizing it via the XXC R-matrix. 


For simplicity, let us hereafter focus on the spin-1 case again. Then the entire Hilbert space is spanned by tensor product of local spaces spanned by the three fundamental vectors |0 ⟩, |1 ⟩, and |2 ⟩. If we set the subset A and B as A = { 0,2 } and B = { 1 }, the XXC Hamiltonian does not change the configuration of the local states |0 ⟩ and |2 ⟩  

    h_XXC   = coshη( |00 ⟩⟨ 00| + |22 ⟩⟨ 22| + |02 ⟩⟨ 02| + |20 ⟩⟨ 20| + |11 ⟩⟨ 11| ) 
       + ( |01 ⟩⟨ 10| + |21 ⟩⟨ 12| + |10 ⟩⟨ 01| + |12 ⟩⟨ 21| ).







  §.§.§ Partially integrable case
 

Now we consider perturbations that break entire integrability but keep integrability in a subspace specified by a given IS. 








The idea is to find integrability-breaking perturbations in such a way that are irrelevant (or vanish) in a given subspace. 

By focusing on systems with nearest-neighbor interactions, we provide two examples of perturbations that violate integrability in the entire Hilbert space but keep integrability in a given subspace. The first one is a perturbation which keeps integrability in the subspace specified by the polarized IS. 






In this subspace, any perturbation in a form of


    h^ pol(α_ d1,α_ d2,α_ o1,α_ o2,ζ)    = α_ d1 |02 ⟩⟨ 02| + α_ d2 |20 ⟩⟨ 20| + α_ o1 |02 ⟩⟨ 20| + α_ o2 |20 ⟩⟨ 02| 
        
          + 2coshζ_1 ∑_s ∈{0,2} |ss ⟩⟨ ss|
        + 2coshζ_2 ∑_s ∈{0,2}( |s1 ⟩⟨ s1| + |1s ⟩⟨ 1s| )


does not violate integrability. Note that these perturbations also preserve the spin-flip invariance. It is easy to check that the first four interactions are irrelevant in the subspace of the polarized IS, since they vanish unless the configuration would include adjacent 0 and 2, which never appears in the polarized subspace. 
On the other hand, the interactions in the second line of Eq. (<ref>) act as a uniform external magnetic field in the projected space P_ polℋ∖{|0⟩,|2⟩}^N. 


Thus, the entire Hamiltonian acts in the projected space as the spin-1/2 XXZ Hamiltonian with the shifted anisotropy coshη→coshη + coshζ_1 - coshζ_2 and the uniform external magnetic field,








    h_j,j+1^XXC + h_j,j+1^ polP_ polℋ∖{|0⟩,|2⟩}^N⟼    σ_j^+ σ_j+1^- + σ_j^- σ_j+1^+ + 1/2(coshη + coshζ_1 - coshζ_2)  (σ_j^z σ_j+1^z) 
       + 1/2coshζ_1  (σ_j^z  + σ_j+1^z) 
        + 1/2 (coshη + coshζ_1 + coshζ_2),


which is integrable. 

The second example is a perturbation which keeps integrability in the subspace of alternating IS.







In this subspace, any of the following perturbations does not violate integrability,


    h^ alt(β_ d1,β_ d2,β_ o1,β_ o2,ζ)    = β_ d1 |00 ⟩⟨ 00| + β_ d2 |22 ⟩⟨ 22| 
        + β_ o1 |00 ⟩⟨ 22| + β_ o2 |22 ⟩⟨ 00| 
          + 2 coshζ_1 (|02 ⟩⟨ 02| + |20 ⟩⟨ 20|)
        + 2coshζ_2 ∑_s ∈{0,2}( |s1 ⟩⟨ s1| + |1s ⟩⟨ 1s| ).


One can see that the first four terms are irrelevant, since they always vanish unless the state includes adjacent 0s or 2s, which never show up in the alternating subspace. On the other hand, the last two terms act as a uniform external magnetic field 
in the projected space P_ altℋ∖{|0⟩,|2⟩}^N. 

Thus, the entire Hamiltonian acts in this projected space as the spin-1/2 XXZ model with the shifted anisotropy coshη→coshη + coshζ_1 - coshζ_2 and the uniform external magnetic field,








    h_j,j+1^ XXC + h_j,j+1^ altP_ altℋ∖{|0⟩,|2⟩}^N⟼    σ_j^+ σ_j+1^- + σ_j^- σ_j+1^+ + 1/2(coshη + coshζ_1 - coshζ_2)  (σ_j^z  σ_j+1^z) 
       + 1/2coshζ  (σ_j^z  + σ_j+1^z) 
        + 1/2 (coshη + coshζ_1 + coshζ_2),


which is again integrable. 


0
By concentrating on the nearest neighbor interactions, the perturbation 

    h'_ alt   = h^00_00 |00 ⟩⟨ 00| + h^22_22 |22 ⟩⟨ 22| 
        + h^02_02 |02 ⟩⟨ 02| + h^20_20 |20 ⟩⟨ 20| 
    
           + h^22_00 |00 ⟩⟨ 22| + h^00_22 |22 ⟩⟨ 00|

keeps integrability in the subspace of the alternating irreducible strings … 0 2 0 2 0 2 …. 

In this subspace, the first two and the last two terms always act as zero, and therefore, the system does not feel the integrability-breaking perturbation in the subspace of the alternating irreducible string. That is, the Hamiltonian after perturbed with these four terms is effectively the XXC model in this subspace. 
On the other hand, the Hamiltonian does not act as the XXC model, after adding the third and fourth terms, even in this subspace. However, we can observe that integrability is still there by “identifying" the states |0 ⟩ and |2 ⟩. Indeed, the Hamiltonian with all terms in (<ref>) is regarded as the spin-1/2 XXZ model after identifying the two states |0 ⟩ and |2 ⟩. We checked this survived integrability as emergence of the embedded spin-1/2 XXZ spectrum in the full spectrum of the perturbed XXC model (see Appendix <ref>). 


Note that the third and fourth terms in (<ref>) violate not only integrability but also the spin-flip invariance, and consequently, the HSF structure of the entire Hilbert space. With these perturbations, total magnetization is no longer a conserved quantity. It is also worth notifying that both of the models (<ref>) and (<ref>) coincide with the t-J_z model <cit.>, a canonical model that exhibits the HSF. 

Another remarkable fact is that partial solvability we discussed above is robust against site-dependent perturbations. For example, one can keep the system integrable in the subspace specified by the alternating IS even when one adds different perturbations on different sites, as long as the perturbations are written in the form of Eq. (<ref>). We will come back to this point later in the discussion of partial solvability for open quantum systems. 





Although we have focused on nearest-neighbor interactions so far, 

there exist longer-range interactions that keep integrability in the subspace specified by an IS.

For instance, the following three-body interactions 
do not violate integrability in the subspace specified by the period-three triplet IS (… 0 0 2 0 0 2 …):

    h'_ tri = h^000_000 |000 ⟩⟨ 000| + h^*22_*221⊗ |22 ⟩⟨ 22| + h^22*_22* |22 ⟩⟨ 22| ⊗1 + h^202_202 |202 ⟩⟨ 202|.




In this way, integrability in the subspace specified by the IS with period-p seems to hold under a certain choice of p-body interactions. Also, the projector onto the subspace encoded by the period-three triplet IS can be written in an (almost) matrix product form with a bond dimension three,

    P_ tri 
        = ∑_0pt(α_1,β_1),…,(α_N,β_N)∈{0,1,2}^2 tr_a (M_α_1,β_1⋯ M_α_N,β_N)
        |α_1 …α_N ⟩⟨β_1 …β_N|
        - 2 ⊗_j=1^N (|1 ⟩⟨ 1|)_j, 
       M_0,0 = S_a^+,    M_1,1 = 1_a,    M_2,2 = (S_a^-)^2,    M_α,β = 0    (α≠β),

where S_a^± and 1_a represent the spin-1 operators and the three-by-three unit matrix acting non-trivially on the auxiliary space, respectively.







  §.§.§ Matrix product operator symmetry
 


In the previous subsection, we have observed that the integrable spin-1/2 XXZ model can be embedded in one of the fragmented Hilbert space of the XXC model specified by a certain IS. Therefore, it is obviously possible to apply the Bethe ansatz method to construct partially conserved quantities, which are conserved only in the integrable subspace but not in the entire Hilbert space. On the other hand, it has been proposed that a partially solvable model is characterized by “the matrix product operator (MPO) symmetry" <cit.>, which implies the existence of conserved quantities in the matrix product forms with fixed bond dimensions. Surprisingly, these are the conserved quantities not only in the integrable subspace but also in the entire Hilbert space. Since they have the matrix product forms with finite bond dimensions, they exhibit small entanglement entropies. 
Concrete examples have been displayed for the XXC model <cit.> under the assumption that the conserved quantity is expressed by the MPO form,

    T =  tr_a (L_a,N… L_a,1), 
       L_a,n = 1⊗ |1 ⟩⟨ 1| + ∑_s,t=0,2 L^(s,t)⊗ |s ⟩⟨ t|,

in which 1 and L^(s,t) are two-by-two matrices acting on the auxiliary space. 

The key relation in proving that the MPO (<ref>) commutes with the Hamiltonian is the local divergence relation,


    [h_j,j+1, L_a,j+1 L_a,j] = M_a,j+1 L_a,j - L_a,j+1 M_a,j,

where M_a,n is another two-by-two matrix. 
0
The mechanism that the MPO (<ref>) commutes with the Hamiltonian comes from the local divergence relation: 

    [h_x,x+1, L_a,x+1 L_a,x] = M_a,x+1 L_a,x - L_a,x+1 M_a,x,

where M is another operator which acts on the tensor product of the two-dimensional auxiliary space and the three-dimensional local physical space. Especially for the Hamiltonian (<ref>) and the MPOs (<ref>) and (<ref>), the local divergence relation (<ref>) is satisfied with M = 0. 


For the XXC model with the nearest-neighbor perturbations (<ref>), only M = 0 solves the divergence relation, which is included in the class of “the commutant algebra" discussed in <cit.>. 
Two kinds of L can be found as the solution to (<ref>): The first one is the diagonal MPO,

    L^(0,0) = [ x y; y z ],   
        L^(2,2) = [ u 0; 0 v ],

where x,y,z and u,v are free parameters. The second one is the non-diagonal MPO, 

    L^(0,2) = (L^(2,0))^† = γσ^-, 
       L^(0,0) = [ α 0; 0 δ ],   
        L^(2,2) = [ β 0; 0 ε ],

where α, β, γ, δ, ε are free parameters. 


Note that the MPO symmetry with M ≠ 0 has been found for the XXC model with longer-range interactions <cit.>, which is not included in the class of the commutant algebra <cit.>. 






§ PARTIALLY SOLVABLE OPEN SPIN CHAINS



Our main focus in this paper is partially solvable quantum systems coupled to boundary dissipators, in which the steady state and some eigenmodes are again solvable. In this section, we show two mechanisms to construct those models: 
The first one is an rSGA-induced partially solvable system which remains to be partially solvable even in the presence of boundary dissipators under a certain condition. The second one is the HSF-induced partially solvable system, whose HSF structure is partially inherited to the boundary dissipative system. 
These are thus examples of partially solvable boundary dissipative systems induced by partial solvability of the Hamiltonian, whose solvable states are robust against the boundary dissipators. 

In addition to robustness of partial solvability against boundary dissipators, there are several important questions including what are characteristic features of partially solvable eigemodes, how relaxation processes in the solvable subspace differ from those in a generic case, and how the partially solvable eigenmodes are experimentally realizable. 


In order to discuss these points, we consider a partially solvable spin-1 chain coupled to boundary dissipators, whose density matrix ρ evolves according to the Liouvillian ℒ in the GKSL equation,

    d/dtρ(t) = ℒ(ρ) 
        = -i [H, ρ] + ∑_αγ_α𝒟_α(ρ), 
       𝒟_α(ρ) = A_αρ A_α^† - 1/2{A_α^† A_α, ρ}.

Here H is the system's Hamiltonian and A_α is a quantum jump operator acting on the physical space. We assume that the jump operator A_α non-trivially acts only on the first and Nth sites, representing the effect of a boundary dissipator with dissipation rates γ_α. 
In the following discussion, we focus on two different Hamiltonians: the AKLT-type Hamiltonian (<ref>) 

and the XXC Hamiltonian (<ref>).




0
In the solvable subspace of the Liouvillian, i.e. equivalently the effective Hamiltonian (<ref>), the dissipation terms always act as zero, and therefore, the dissipators are irrelevant in the solvable subspace. That is, any state in this subspace is a so-called “dark state" <cit.>, which does not feel the dissipators. 
In the following subsections, we propose three mechanisms for emergence of exactly solvable eigenmodes induced by bulk partial solvability. 





 §.§ rSGA induced solvable eigenmodes

The first example in which partial solvability is robust against boundary dissipators is given by an rSGA-induced partially solvable system. Let us consider a system in which a zero-energy state |Ψ_0 ⟩ is analytically known. The rSGA solvability is characterized by the existence of a spectrum generating operator Q^† that satisfies the rSGA (<ref>) with the Hamiltonian in a subspace W spanned by states constructed by applying Q^† to the zero-energy state,

    W = { |Ψ_n ⟩}_n,   
        |Ψ_n ⟩ = (Q^†)^n |Ψ_0 ⟩.

Thus, the states |Ψ_n ⟩ are exactly solvable energy eigenstates of the Hamiltonian with the eigenenergies ℰn. 
It often occurs that the solvable excited states admit the single-mode quasiparticle description,

    Q^† = ∑_x=1^N e^ikx q^†_x,

in which q_x^† is a local operator non-trivially acting on the xth site. 

Based on these, it is natural to ask whether the solvable states can survive even when quasiparticles are injected from both of the edges of the system.
Such a situation is realized by taking the boundary quantum jump operators as 

    A_ L = q_1^†,   
        A_ R = q_N^†.







The density matrix ρ_nn for the solvable energy eigenstate |Ψ_n ⟩, as it commutes with the Hamiltonian by definition (i.e., [ρ_nn, H]=0), becomes the steady state of the GKSL equation (<ref>) if 

    𝒟_α(ρ_nn) = 0,   ∀α, 
       ρ_nn = |Ψ_n ⟩⟨Ψ_n|.


The pure state ρ_nn that satisfies this condition together with the commutativity with the Hamiltonian is known as “the dark state", which has been introduced in the context of atomic physics and optics <cit.>. 

The AKLT-type model (<ref>) is one of the examples whose solvable energy eigenstates satisfy the dark-state condition (<ref>) in the presence of the boundary quasiparticle dissipators.


The zero-energy eigenstate is given by the matrix product state, as has been explained in Sec. <ref>, but with the boundary deformation,

    |Ψ^(v_ L, v_ R)_0 ⟩ 
        = ∑_m_1,…,m_N ∈{0,1,2}_a⟨ v_ L| A_m_1⋯ A_m_N |v_ R⟩_a · |m_1… m_N ⟩,


since the open boundary condition is imposed on the system. The boundary vectors |v_ L,R⟩∈ V_a =  span{ |↑⟩, |↓⟩} in the auxiliary space must be properly chosen in order for (<ref>) to be the zero-energy eigenstate. Especially when the model is frustration-free, as we consider in this paper, there are four degenerate zero-energy eigenstates, as no constraint is imposed on the boundary vectors. 
A tower of the solvable energy eigenstates are then independently constructed on top of each of the four degenerate zero-energy states by applying the spin-2 magnon creation operator Q^† (<ref>). That is, the solvable subspace W under the open boundary condition is composed of four separate subspaces specified by the boundary vectors,

    W = W^(↑,↑)⊕ W^(↑,↓)⊕ W^(↓,↑)⊕ W^(↓,↓), 
       W^(v_ L,v_ R) =  span{ |Ψ_n^(v_ L,v_ R)⟩}_n,    
        |Ψ_n^(v_ L,v_ R)⟩ = (Q^†)^n |Ψ_0^(v_ L,v_ R)⟩.


We set the boundary quasiparticle dissipators in such a way that the quasiparticles are coming into the system from both of the ends,

    A_ L = (S_1^+)^2,    
        A_ R = (S_N^+)^2.







With these dissipators, one of the solvable subspaces W^(↑,↓) satisfies the dark state conditions (<ref>),

    A_α |ψ^(↑,↓)⟩ = 0,   α∈{ L,  R}, 
       |ψ^(↑,↓)⟩∈ W^(↑,↓).

(See Appendix <ref> for the proof.) That is, by denoting the Liouvillian with the AKLT-type Hamiltonian and the dissipators (<ref>) by ℒ_ AKLT, any diagonal density matrix in the subspace W^(↑,↓) 
becomes a steady state of the GKSL equation,

    ℒ_ AKLT(ρ_ diag^(↑,↓)) = 0, 
       ρ_ diag^(↑,↓) = ∑_n p_n |Ψ_n^(↑,↓)⟩⟨Ψ_n^(↑,↓)|,   ∑_n p_n = 1,   p_n ≥ 0, ∀ n.


Note that, in the case of the perturbed spin-1 XY model, which is another model with the hidden rSGA, there never exist such a solvable energy eigenstate that satisfies the dark-state condition (<ref>).  

The fact that any state in the subspace W^(↑,↓) is the dark state (<ref>) indicates that any density matrix in W^(↑,↓)⊗ (W^(↑,↓))^†, even if it contains off-diagonal elements, becomes an eigenmode of the GKSL equation. Suppose that we have an off-diagonal element ρ_m,n^(↑,↓) = |Ψ_m^(↑,↓)⟩⟨Ψ_n^(↑,↓)| in a given density matrix. This is a dark state from the previous statement, and moreover, its commutator with the Hamiltonian is proportional to itself,

    [H, ρ_m,n^(↑,↓)] = 2(m-n) h^00_00ρ_m,n^(↑,↓),

which indicates that ρ_m,n^(↑,↓) is an eigenmode of the GKSL equation. The relation (<ref>) together with the dark-state condition (<ref>) leads to the restricted spectrum generating algebra for the Liouvillian in the subspace W^(↑,↓)⊗ (W^(↑,↓))^†⊂ℋ⊗ℋ^†,

    ℒ_ AKLT(ρ^(↑,↓)_m,n)
        = -2i (m-n) h^00_00ρ^(↑,↓)_m,n,

giving an equally-spaced spectrum along the imaginary axis with the interval 2h^00_00 embedded in the full spectrum of ℒ_ AKLT. 









The hidden rSGA structure (<ref>) of the Liouvillian evokes us the persistent oscillations observed for rSGA-induced partially solvable isolated quantum systems <cit.>. Indeed, if we choose an initial state in the subspace W^(↑,↓),

    |ψ(0) ⟩ = ∑_n a_n |Ψ_n^(↑,↓)⟩∈ W^(↑,↓),

it is easy to show that persistent oscillations are observed for an observable O in a long-time scale,

    ⟨ O(t) ⟩∼∑_n ≤ m 2cos(2(m-n)h^00_00 t) a_m a_n  Re  O_nm   (t→∞).

This indicates that the system prepared in the solvable subspace never relaxes to any steady state. 
Even if the initial state is generic, the long-lived oscillations survive as long as the initial state has a large enough overlap with the subspace W^(↑,↓). 

0
mechanism we propose in this paper is the partially solvable Liouvillians associated with the rSGA. The mechanism of partial solvability for this kind of models is explained by the robust rSGA under the boundary dissipators. That is, the bulk Hamiltonian has partial solvability induced by the rSGA, and the rSGA structure is not violated by the boundary dissipators. 
This situation is realized if the solvable subspace W is a subspace of the kernel of all the dissipation operators 𝒟_α. That is, the following relations must hold for all |ψ⟩∈ W

    H |ψ⟩ = E |ψ⟩,     𝒟_α(|ψ⟩⟨ψ|) = 0,  ∀α.

Similar steady states have been known by the name “dark states", first introduced in <cit.> and discussed later in many examples <cit.>, which are defined by slightly stronger conditions

    H |ψ⟩ = E |ψ⟩,     
        A_α |ψ⟩ = 0,  ∀α.

For this reason, we call the state which satisfies (<ref>) “the generalized dark state". 

Several examples of the dark states have been found for the systems in which every site is coupled to a single dephasing operator <cit.>. These models exhibits the persistently oscillating eigenmodes, which are the characteristic behavior of the models equipped with the rSGA. 

On the other hand, the system which is coupled to dissipators only at the edges can admit the generalized dark states. Such generalized dark states are obtained for the partially solvable spin chain with frustration-free interactions in the spin-2 magnon baths at the edges. The examples include the AKLT model and its generalizations (<ref>). 
These models exhibit spin-2 magnon excitation (<ref>) on the top of the zero-energy matrix product state (<ref>), which form the exactly solvable invariant subspace W of the Hamiltonian (<ref>). 

Suppose that the partially solvable Hamiltonian with the rSGA is coupled to boundary dissipators given by (S_1^+)^2 and (S_N^+)^2. 
The first relation for the generalized dark state condition (<ref>) is satisfied if we choose |ψ⟩ = (Q^†)^n |Ψ_0 ⟩ as it is the energy eigenstates, while the second relation requires 

    (S_1^+)^2 (Q^†)^n |Ψ_0 ⟩⟨Ψ_0| Q^n (S_1^-)^2 - 1/2{ (S_1^-)^2 (S_1^+)^2,  (Q^†)^n |Ψ_0 ⟩⟨Ψ_0| Q^n } = 0, 
       (S_N^+)^2 (Q^†)^n |Ψ_0 ⟩⟨Ψ_0| Q^n (S_N^-)^2 - 1/2{ (S_N^-)^2 (S_N^+)^2,  (Q^†)^n |Ψ_0 ⟩⟨Ψ_0| Q^n } = 0.

These are satisfied by one of the choices for the boundary vectors ** which produces four degenerate zero-energy eigenstates (<ref>). 

It is important for obtaining the generalized dark states that the vacuum of quasiparticles are the non-trivial states. In fact, the solvable subspace of the perturbed spin-1 XY model consists of the vacuum, which is the fully-polarized state (<ref>), and the spin-2 magnon excitation on the top of the vacuum (<ref>), and this subspace is violated by the quasiparticle baths at the edges. 
The non-trivial structure of the zero-energy state for the AKLT model is understood by the deformed symmetry structure <cit.>, which connects the non-trivial zero-energy state with the fully-polarized state. Thus, it is necessary for the solvable subspace of the bulk Hamiltonian is associated with the deformed symmetry structure that to be robust against the quasiparticle baths at the edges. 

As the dissipators are irrelevant in the solvable subspace of the rSGA-induced partially solvable models, persistent oscillating modes exist for the Liouvillian, as a consequence of the equally-spaced spectrum of the bulk Hamiltonian. These oscillating modes are characterized by 

    [H, ρ_m,n] = -ω (m-n) ρ_m,n, 
       𝒟_α(ρ_m,n) = 0,   ∀α,

leading to 

    ρ(t)    = ∑_n c_n^2 |Ψ_n ⟩⟨Ψ_n| 
       + ∑_n ≠ m c_n c_m { e^iω (m-n) t |Ψ_n ⟩⟨Ψ_m| + e^-iω (m-n) t |Ψ_m ⟩⟨Ψ_n| },

and thus to the persistent oscillation of local observables, if the initial state has a decent overlap with the quasiparticle excitation states, indicating that the system never reaches the relaxation state. 

Similar phenomena have been observed for the system with the dephasing operators coupled to every site <cit.>. Here, we observed that the non-relaxing states exist even when the system is coupled to the boundary dissipators, each of which dopes quasiparticles. 





 §.§ Numerical simulation


To see how the presence of the rSGA-induced solvable eigenmodes affects observables in partially solvable open spin chains,
here we perform numerical simulations for the GKSL equation (<ref>).
We take the generalized AKLT model H_ AKLT=∑_j=1^N-1 h_j,j+1^ AKLT (<ref>) as the bulk Hamiltonian
with a finite number of lattice sites N. We specifically focus on the case of the original AKLT model, i.e., by choosing
h_11^11/h_00^00=2/3, a_0=-√(2)a_1=-a_2=√(2/3).
The dissipators are given by spin-2 creation operators (<ref>)
acting on each end of the chain. The explicit form of the GKSL equation that we solve in this section is:

    d/dtρ   =
    -i[H_ AKLT,ρ]+∑_α= L,Rγ_α( A_αρ A_α^† - 1/2{A_α^† A_α, ρ}),
    
    
    H_ AKLT   =
    ∑_j=1^N-1[
     S_j · S_j+1+1/3 ( S_j · S_j+1)^2
    ],
    
    
    A_ L   =
    (S_1^+)^2,
      
    A_ R
    =
    (S_N^+)^2,

which is numerically simulated by the quantum trajectory method <cit.> together with the exact diagonalization.

The initial state is chosen to be a product state of the spin-1 chain, whose spin configuration is nearly a Néel state |N, S^z(0)⟩,
depending on the values of N and the initial total S^z(0):

    |0 2 0 2 ⋯ 0 2 0 2 ⟩ : 
         
    
    |0 2 0 2 ⋯ 0 2 0 1 ⟩ : 
         
    
    |0 2 0 2 ⋯ 0 2 0 2 1 ⟩ : 
         
    
    |0 2 0 2 ⋯ 0 2 0 2 0 ⟩ :

We identify the state labels 0, 1, 2 with the spin configuration ↑, 0, ↓, respectively. Among the initial states |N, S^z(0)⟩ considered here, those with S^z(0)=1 have an overlap with the states in the subspace W^(↑,↓), since they have ↑ spins at both ends of the chain after removing S_j^z=0 spins. Hence, we expect that the rSGA-induced solvable eigenmodes appear in those cases with S^z(0)=1.

In Fig. <ref>, we show the time evolution of the local magnetization ⟨ S_j^z⟩ for several initial conditions
((a) N=8, S^z=0, (b) N=8, S^z=1, (c) N=9, S^z=0, (d) N=9, S^z=1)
with γ_ L=γ_ R=1. In the cases of S^z(0)=0 [Fig. <ref>(a), (c)], the local magnetization gradually reaches 
the steady state value without oscillations, while in the cases of S^z(0)=1 [Fig. <ref>(b), (d)]
the local magnetization clearly shows long-lived coherent oscillations. In all the cases above, the magnetization does not approach the maximum value, meaning that the steady state does not correspond to the trivial all up states (i.e., |00⋯ 0⟩=|↑↑⋯↑⟩). We can also see that the magnetization at the boundary j=1, N takes a relatively larger steady-state value as compared to the bulk part, which may be due to the effect of the spin injection at the boundary.



In order to understand the role of the solvable eigenmodes in the GKSL equation, we plot the ratio of the number of trajectories for each S^z measured by ⟨ P_S^z⟩ in Fig. <ref>, where P_S^z is the projection operator onto the corresponding subspace with fixed S^z. The parameters are the same as in Fig. <ref>. When S^z(0)=0 [Fig. <ref>(a), (c)], the number of trajectories having S^z=0 quickly decays to zero, while those with S^z≠ 0 grow subsequently (those with lower S^z grow faster). Since S^z can change by 2 due to the spin-2 injection at the boundaries, S^z only takes values of even integers (that should not exceed the system size N). In the long-time limit, 
the trajectories with S^z=6 and 8 survive, and the others seem to vanish.
This is in sharp contrast to the cases for S^z(0)=1 [Fig. <ref>(b), (d)],
where the number of trajectories with arbitrary S^z can survive in the long-time limit. This is in consistent with the fact that there is a tower of dark states with S^z=1,3,5,… (as shown in Eq. (<ref>)), which can be accessed from the initial state with S^z(0)=1 having an overlap with the states in the subspace W^(↑,↓). Hence the steady state remains to be far from the trivial all up states (|00⋯ 0⟩=|↑↑⋯↑⟩) even in the presence of the spin-2 injection. The steady states for S^z(0)=1 are also different from the all up state, since there is an additional dark state with S^z=S_ max^z-2 (S_ max^z is the maximum S^z in a finite chain with length N). In the thermodynamic limit, this dark state will become indistinguishable from the all up state.
In this way, there is a clear difference between the dynamics starting from S^z(0)=0 and S^z(0)=1 rooted in the presence of the solvable eigenmodes in the GKSL equation.



The long-lived coherent oscillations observed in Fig. <ref>(b), (d) are not due to the solvable dark states, since each trajectory has a single value of S^z at each time, which allows for realization of a single dark state at each time in each trajectory. One cannot have quantum mechanical superposition of different dark states,
which might cause coherent oscillations due to interference among dark states. To find the origin of the oscillations, we plot the imaginary part of the eigenvalues of the non-Hermitian Hamiltonian H^ eff_ AKLT=H_ AKLT-i/2∑_α= L,Rγ_α A_α^† A_α 
with S^z=0 in Fig. <ref>(a) and S^z=1 in Fig. <ref>(b). In the former case (S^z=0), there is no eigenstate having an eigenvalue with zero imaginary part.
All the eigenstates have nonzero imaginary parts, which are forming continuous spectra. 
In the latter case (S^z=1), on the other hand, there is one eigenstate having an eigenvalues with zero imaginary part, which is not shown in Fig. <ref>(b) where we take a log scale in the vertical axis. This corresponds to one of the solvable eigenmodes in the GKSL equation. On top of that, we find another eigenstate having an eigenvalue with a nonzero but very small imaginary part, which is separated from those of the other eigenstates. This eigenstate has the second smallest real part of the eigenvalues. Although it is not exactly a dark state, it may create oscillations with lifetime longer than the maximum time in Fig. <ref>. In fact, we confirm that this eigenstate has a strong overlap with the initial state close to the Néel state for S^z(0)=1. If one waits for sufficiently long time, we expect that the coherent oscillations will vanish eventually.



We also performed numerical simulations for other choices of the model parameters away from the solvable regions. The results show similar behaviors (e.g., the number of trajectories with S^z<S_ max^z-2 decay to zero) as in the case of S^z(0)=0 for the AKLT model with the boundary dissipators. Again we confirm the role of the rSGA-induced solvable eigenmodes in the GKSL equation, which support nontrivial steady states in the presence of dissipations.





 §.§ HSF-induced solvable eigenmodes


The second example of robust partial solvability against boundary dissipators can be found for systems with the HSF. 


The HSF of the Liouvillian has been discussed in the context of the commutant algebra <cit.>. 

In this subsection, on the other hand, we discuss partial solvability induced by the HSF structure of the Hamiltonian for boundary dissipative systems, which have not been considered before. 



In order to discuss the HSF for Liouvillians, it is useful to work on the TFD vector expression <cit.>, which is realized by the isomorphism,

    φ:    ρ = ∑_m,n p_m,n |m ⟩⟨ n| 
        ↦ |ρ⟩⟩ = ∑_m,n p_m,n |m ⟩⊗ |n ⟩^*.

In the TFD expression, the Liouvillian is expressed as a non-Hermitian Hamiltonian acting on the doubled Hilbert space ℋ⊗ℋ^*,

    d/dt |ρ(t) ⟩⟩ = -i H |ρ(t) ⟩⟩, 
       H = H ⊗1 - 1⊗^tH + i ∑_αγ_α( (A_α⊗ A_α^*) - 1/2 (A_α^† A_α⊗1 + 1⊗^tA_α A_α^*) ).

For boundary dissipative systems, the quantum jump operators A_α non-trivially act only on the first and/or the Nth site. In this subsection, we consider the Hamiltonian with the HSF that is solvable at least in one of the fragmented subspaces. 

0
One may notice that this non-Hermitian Hamiltonian consists of the Hermitian part given by the first two terms and the non-Hermitian part given by the rests. Thus, a state |ψ⟩⟩∈ V ⊗^tV can be the zero-energy eigenvector of H, by proving the steady state of ℒ, if |ψ⟩⟩ separately vanishes by the first and second parts

    (H ⊗1 - 1⊗^tH) |ψ⟩⟩ = 0, 
       ( (A_α⊗ A_α^*) - 1/2 (A_α^† A_α⊗1 + 1⊗^tA_α A_α^*) ) |ψ⟩⟩ = 0,   ∀α.

These conditions are sufficiently satisfied if the state |ψ⟩⟩ is a so-called “dark state" i.e. an energy eigenstate and killed by all quantum jump operators

    H |ψ⟩ = E |ψ⟩, 
       𝒟_α(|ψ⟩⟨ψ|) = 0,   ∀α.





Throughout this subsection, we consider the (perturbed) XXC Hamiltonian (<ref>) that exhibits both the HSF and (partial) integrability. 








Suppose that two kinds of dissipators are coupled to each end of the XXC spin chain,

    A_ L,+ = (S_1^+)^2,    A_ L,- = (S_1^-)^2,    
        A_ R,+ = (S_N^+)^2,    A_ R,- = (S_N^-)^2,

with the coupling strengths controlled by the dissipation rates γ_ L,+, γ_ L,-, γ_ R,+, and γ_ R,-. 

With these dissipators, the Liouviilian for the boundary dissipative XXC spin chain ℒ_XXC is effectively written as the spin chain having twice the length 2N than the original chain (Fig. <ref>),

    H_XXC = ∑_j=1^N-1 h^(XXC)_j,j+1 + h^( b,R)_N,N+1 
        - ∑_j=N+1^2N-1 h^(XXC)_j,j+1 + h^( b,L)_1,2N.


Here we have used the transpose invariance and the inversion symmetry of the XXC Hamiltonian (<ref>),

    ^t(∑_j=1^N-1 h_j,j+1^XXC) = ∑_j=1^N-1^th_j,j+1^XXC = ∑_j=1^N-1 h_j,j+1^XXC, 
       U_ I(∑_j=1^N-1 h_j,j+1^XXC) U_ I 
        = ∑_j=1^N-1 U_ I h_j,j+1^XXC U_ I
        = ∑_j=1^N-1 h_N-j+1,N-j^XXC = ∑_j=1^N-1 h_j,j+1^XXC,

in which the operator U_ I reflects the spin chain with respect to its center,

    U_ I: ℋ→ℋ,   |s_1,s_2,…,s_N ⟩↦ |s_N,…,s_2,s_1 ⟩.



The effects of the boundary dissipators can be written in terms of interactions between the first and 2Nth sites and between the Nth and (N+1)th sites,

    h^( b,α)   = iγ_α,+( |00 ⟩⟨ 22| - 1/2 (|2 ⟩⟨ 2| ⊗1 + 1⊗ |2 ⟩⟨ 2|) ) 
          + iγ_α,-( |22 ⟩⟨ 00| - 1/2 (|0 ⟩⟨ 0| ⊗1 + 1⊗ |0 ⟩⟨ 0|) ),   α∈{ R,L},

each of which represents incoming and outgoing quasiparticles.

When all the dissipation rates are set to zero, the effective Hamiltonian (<ref>) simply consists of the two decoupled XXC Hamiltonians. Therefore, the HSF structure is obtained in the doubled Hilbert space according to the IS consisting of configurations of 0s and 2s. Moreover, this doubled XXC Hamiltonian is obviously integrable, since each of the XXC spin chains is integrable. 

On the other hand, both the HSF structure and entire integrability are broken by the presence of boundary dissipation terms. However, one may notice that the subspace specified by the alternating IS survives as an invariant subspace of the entire effective Hamiltonian (<ref>) even in the presence of the boundary dissipation terms, since the nearest-neighbor interactions in the boundary dissipation terms (<ref>), i.e. the IS violating terms, are irrelevant in this subspace. 

In the original Liouvillian expression, this means that any state with the alternating IS never goes out from the other subspaces, nor any state in the other subspaces never come into the subspace with the alternating IS.









  §.§.§ Integrable subspaces

The diagonal terms in (<ref>) can be regarded as boundary magnetic fields on each of the two XXC Hamiltonians, and therefore the subspace specified by the alternating IS is the integrable subspace of the effective Hamiltonian (<ref>) if these diagonal terms provide the integrable boundaries of the XXC model. 





The integrable boundary conditions for the XXC model have been discussed in <cit.>. 
When boundary magnetic fields are imposed only in the z-direction, the Hamiltonian is integrable if the boundary terms are given by 

    H_ bXXC   = ∑_j=1^N-1h^XXC_j,j+1 
        + sinhηξ_- · (S_1^z)^2
        - sinhηξ_+ · (S_N^z)^2 
          + 1/2sinhη (ξ_+ - ξ_-).

Here ξ_- and ξ_- are arbitrary complex parameters. 




0
The R-matrix for the XXC model is given in (<ref>). The newly introduced matrix K representing the reflection at the boundaries  contains the parameter ξ, which determines strength of the diagonal boundary magnetic fields. The R-matrix and K-matrix for the XXC model can be found in <cit.>. 

With these observations, the effective Hamiltonian (<ref>) becomes integrable in the subspace specified by the alternating IS when the quasiparticle incoming and outgoing rates are the same at each end,

    γ_α,+ = γ_α,-,   α∈{ L,  R}.


Although we set H as the integrable XXC Hamiltonian so far, it is also possible to replace it with the perturbed XXC Hamiltonian by adding site-dependent perturbations on the bulk without violating partial solvability of the Liouvillian in the subspace with the alternating IS. Such a perturbation (<ref>) modifies the bulk interactions in the effective Hamiltonian (<ref>) as 

    h^XXC_j,j+1→ h^XXC_j,j+1 + h_j,j+1^ alt(β^(j)_ d1,β^(j)_ d2,β^(j)_ o1,β^(j)_ o2,ζ^(j)),    j=1,…,N-1, 
       h^XXC_j,j+1→ h^XXC_j,j+1 + h_j,j+1^ alt(β^(j)_ d1,β^(j)_ d2,β^(j)_ o2,β^(j)_ o1,ζ^(j)),    j=N+1,…,2N-1,

both of which are apparently in the form of (<ref>), implying that the effective Hamiltonian stays partially solvable in the subspace with the alternating IS after introducing these bulk perturbations. 








0
We found four kinds of subspaces in each of which the effective Hamiltonian (<ref>) with the dissipators (<ref>) becomes different integrable Hamiltonians. The projectors onto these subspaces are essentially the matrix product projector onto the subspace of the alternating irreducible strings (<ref>), but have different boundary conditions for different subspaces. 

Let us introduce two different kinds of boundary projectors which non-trivially act on the first or last site

    B_0, L = (|0 ⟩⟨ 0| + |1 ⟩⟨ 1|) (⊗1)^N-1,   
        B_2, L = (|2 ⟩⟨ 2| + |1 ⟩⟨ 1|) (⊗1)^N-1, 
       B_0, R =  (1⊗ )^N-1 (|0 ⟩⟨ 0| + |1 ⟩⟨ 1|),   
        B_2, R =  (1⊗ )^N-1 (|2 ⟩⟨ 2| + |1 ⟩⟨ 1|).

We write the bulk part of the matrix product projector (<ref>) as

    𝒫_ alt = ⊗_ phys( σ_ aux^+ ⊗ |0 ⟩⟨ 0| + 1_ aux⊗ |1 ⟩⟨ 1| + σ_ aux^- ⊗ |2 ⟩⟨ 2| ),

which is related to the projector P_ alt in (<ref>) by P_ alt =  tr_ aux 𝒫_ alt. 
Then, the effective Hamiltonian (<ref>) under any of following projectors 

    P_02,02 =  tr_ aux (B_0,  L𝒫_ alt B_2,  R) ⊗ tr_ aux (B_0,  L𝒫_ alt B_2,  R), 
       P_20,20 =  tr_ aux (B_2,  L𝒫_ alt B_0,  R) ⊗ tr_ aux (B_2,  L𝒫_ alt B_0,  R), 
       P_00,22 =  tr_ aux (B_0,  L𝒫_ alt B_0,  R) ⊗ tr_ aux (B_2,  L𝒫_ alt B_2,  R), 
       P_22,00 =  tr_ aux (B_2,  L𝒫_ alt B_2,  R) ⊗ tr_ aux (B_0,  L𝒫_ alt B_0,  R)

is the two decoupled integrable XXZ model with imaginary boundary magnetic fields. Each projector realizes the different boundary magnetic fields. 

0
written in the matrix product expressions

    𝒫_α_1,α_N= P_α_1,α_N P_α_N,α_1, 
       P_α_1,α_N = ∑_α_1,…,α_N ∈{0,1,2} tr[ A_α_1 A_α_2… A_α_N-1 A_α_N] |α_1 α_2 …α_N-1α_N ⟩⟨α_1 α_2 …α_N-1α_N| 
       - ∑_α_1, α_N ∈{0,2}∑_α_2,…,α_N-1∈{0,1,2} tr[ A_α_1 A_α_2… A_α_N-1 A_α_N] |α_1 α_2 …α_N-1α_N ⟩⟨α_1 α_2 …α_N-1α_N|,

where we used the notations 0 = 2 and 2 = 0. With these projectors, the effective Hamiltonian is mapped to the spin-1/2 XXZ chain under different boundary conditions depending on the choices of α_ L and α_ R. 

by choosing as γ_ L^+ = γ_ L^(-) and γ_ R^+ = γ_ R^(-), the effective Hamiltonian becomes the periodic XXC chain with two impurities which deform the interaction at the N and (N+1)th sites, and the first and 2Nth sites as 

    h_1, 2N = γ_- { |22 ⟩⟨ 00| 
        - 1/2 (|0 ⟩⟨ 0| ⊗1) - 1/2 (1⊗ |0 ⟩⟨ 0|) }, 
       h_N,N+1 = γ_+ { |00 ⟩⟨ 22|  
        - 1/2 (|2 ⟩⟨ 2| ⊗1) - 1/2 (1⊗ |2 ⟩⟨ 2|) },

whose first term acts as zero in the subspace of alternating irreducible strings, while the second and third terms are regarded as imaginary boundary magnetic fields of two decoupled spin chain of length N in the same subspace. 

For instance, the projected Hamiltonian by P_02,02 becomes 

    H_P_02,02𝒫_02,02ℋ∖{0,2}^⊗ N⟼ H^(+)_XXZ(γ_-, L,γ_+, R) ⊗1 - 1⊗ H^(-)_XXZ(γ_-, R,γ_+, L), 
       H^(+)_XXZ(γ_-, L,γ_+, R) = - i/4γ_-, Lσ_1^z + ∑_x=1^N-1( σ_x^+ σ_x+1^- + σ_x^- σ_x+1^+ + 1/2coshη σ_x^z σ_x+1^z ) - i/4γ_+, Rσ_N^z 
       + N-1/2coshη - i/4 (γ_-, L + γ_+, R) 
       H^(-)_XXZ(γ_-, R,γ_+, L) = - i/4γ_-, Rσ_N+1^z - ∑_x=N+1^2N-1( σ_x^+ σ_x+1^- + σ_x^- σ_x+1^+ + 1/2coshη σ_x^z σ_x+1^z ) - i/4γ_+, Lσ_2N^z 
       + N-1/2coshη - i/4 (γ_+, L + γ_-, R)

in the reduced space ⊗^N  span{ |A ⟩, |B ⟩}. We provide the projected Hamiltonian by each of (<ref>) in Table <ref>, which is numerically verified for N = 3 and 4. 













0
through the transfer matrix 

    T(λ) = **.

Note that the transfer matrix constructed in this way has commuting property

    [T(λ), T(μ)] = 0,

which produces many conserved quantities. 
The commuting matrices are constructed as the product of the R-matrices 

    T(λ) =

produces the commuting transfer matrices 

    T(λ) =  tr_0 (R_0N(λ) … R_02(λ) R_01(λ)), 
       [T(λ), T(μ)] = 0,

and subsequently, many conserved quantities as well as the expansion of T(λ) with respect to the spectral parameter λ. 
The energy eigenstates are expressed 

For instance,  the subspace encoded by the fully-polarized irreducible strings … 0000 …, 







  §.§.§ Solvable eigenmodes


The XXC Hamiltonian with the integrable boundaries can be diagonalized via the Bethe ansatz method due to the existence of the R- and K-matrices which solve the Yang-Baxter equation and the reflection relation, respectively. For the perturbed XXC model that becomes integrable only in the subspace specified by a certain IS,

the spectrum and eigenvectors in this subspace can be found by mapping the Hamiltonian restricted in the solvable subspace to the spin-1/2 XXZ model. 

By employing the same strategy explained in Sec. <ref>, the effective Hamiltonian H_XXC (<ref>) in the solvable subspace, namely the subspace specified by the alternating IS, is first mapped to two decoupled XXC spin chains (Fig. <ref>). Then these XXC spin chains can be mapped to the spin-1/2 XXZ chains with the diagonal boundary magnetic fields by identifying the states |0 ⟩ and |2 ⟩,

    H_XXCP_ altℋ∖{ |0⟩,|2⟩}^N⟼ H^(+)_XXZ(γ_ L,γ_ R) ⊗1 - 1⊗ H^(-)_XXZ(γ_ R,γ_ L), 
       H^(+)_XXZ(γ_ L,γ_ R) = - i/4γ_ Lσ_1^z + ∑_j=1^N-1( σ_j^+ σ_j+1^- + σ_j^- σ_j+1^+ + 1/2coshη σ_j^z σ_j+1^z ) - i/4γ_ Rσ_N^z 
       + N-1/2coshη - i/4 (γ_ L + γ_ R), 
       H^(-)_XXZ(γ_ R,γ_ L) = - i/4γ_ Rσ_N+1^z - ∑_j=N+1^2N-1( σ_j^+ σ_j+1^- + σ_j^- σ_j+1^+ + 1/2coshη σ_j^z σ_j+1^z ) - i/4γ_ Lσ_2N^z 
       + N-1/2coshη - i/4 (γ_ L + γ_ R).

Note that the boundary magnetic fields here are “imaginary magnetic fields" with pure imaginary coefficients. By writing sets of energy eigenvalues for these two spin-1/2 XXZ chains as {E^(+)_n_+}_n_+ and {E^(-)_n_-}_n_-, the set of the summed energy eigenvalues {E^(+)_n_+ + E^(-)_n_-}_n_+,n_- is embedded in the full spectrum of the effective Hamiltonian H_XXC. 
Since the map φ is an isomorphism, the spectrum of the Liouvillian ℒ_XXC matches that of the effective Hamiltonian H_XXC. The energy eigenvalues of the double spin chain, {E^(+)_n_+ + E^(-)_n_-}_n_+,n_-, then agree with the eigenvalues of the Liouvillian ℒ_XXC restricted in the solvable subspace. 

The corresponding eigenvectors |E^(+)_n_+⟩⊗ |E^(-)_n_-⟩ thus provide the eigenvectors of the effective non-Hermitian Hamiltonian H_XXC by a map 


    |σ_1 σ_2 …σ_N ⟩∈ℂ^2N↦ |τ_1 τ_2 …τ_N ⟩∈ℂ^3N, 
       τ_j = 1 + θ_j^ω_j,

in which θ_j and ω_j are defined by 

    θ_j = -|1 - σ_j|,     ω_j = 1 - (-1)^∑_k=1^j σ_k,

respectively.
Here θ_j is a parameter that determines whether the jth site is in the local state |1 ⟩ or not, while ω_j counts the number of 0s and 2s between the first and jth site. 

The eigenvectors of the effective non-Hermitian Hamiltonian H_XXC are mapped to the eigenmodes of the Liouvillian ℒ_XXC via the inverse map of the isomorphism (<ref>). 
Therefore, solving the eigenvalue problem for the spin-1/2 XXZ chain under the imaginary boundary magnetic fields tells the eigenmodes for the spin-1 XXC model coupled to boundary dissipators. 





The Bethe-ansatz method is well established for the spin-1/2 XXZ chain even in the presence of boundary magnetic fields <cit.>. The Hamiltonian and conserved quantities are constructed from a series expansion of the transfer matrix, which consists of the R-matrix

    R_ij(λ)    = 
    
    
    
    
    
    
        sinh(λ+ η/2) coshη/2·1_ij + cosh(u+η/2) sinhη/2·σ_i^z σ_j^z 
          + sinhη·( σ_i^+ σ_j^- + σ_i^- σ_j^+ ),

and the K-matrix

    K(λ,ξ) = sinhξcoshλ·1 + coshξsinhλ·σ^z,


which solve the Yang-Baxter equation (<ref>) and the reflection relation,

    R_12(λ_1-λ_2) K_1(λ_1) R_12(λ_1+λ_2) K_2(λ_2)
        = K_2(λ_2) R_12(λ_1 + λ_2) K_1(λ_1) R_12(λ_1-λ_2).

The complex parameter ξ determines the strength of the diagonal boundary magnetic fields, which apprears in the spin-1/2 XXZ Hamiltonian as  

    H_XXZ
           = ∑_j=1^N-1( σ_j^+ σ_j+1^- + σ_j^- σ_j+1^+ + 1/2coshη·σ_j^z σ_j+1^z ) 
          +1/2sinhηξ_- ·σ_1^z  - 1/2sinhηξ_+ ·σ_N^z.

Thus, in order to realize the restricted effective Hamiltonian (<ref>),

we need to choose ξ_± to be pure imaginary for η∈ℝ (i.e., in the gapped regime) or 

to be real for pure imaginary η (i.e., in the gapless regime).

The transfer matrix T(λ) for the open spin chain is “a double-row transfer matrix", which consists of two products of the R-matrices,

    T(λ) =  tr_0 ( K_0(-λ - η,ξ_+)  M_0(λ)  K_0(λ,ξ_-) M_0(λ) ), 
       M_0(λ) = R_0N(λ) … R_01(λ), 
       M_0(λ) = R_10(λ) … R_N0(λ) .

These transfer matrices are mutually commuting,

    [T(λ),  T(μ)] = 0,

for any λ, μ∈ℂ. 
Thus, a series expansion of the transfer matrix,

    T(λ) = exp( ∑_r λ^r/r! Q_r ),

provides a large number of conserved quantities Q_r. By cumbersome but straightforward calculations, one can confirm that the XXZ Hamiltonian with the boundary magnetic fields (<ref>) is obtained 

from Q_1 up to a constant <cit.>,

    H_XXZ   = (2 sinhξ_+ sinhξ_- coshη  (sinhη)^2N-1)^-1·d/dλ T(λ) |_λ = 0
        - (sinhη)^2 + N (coshη)^2/coshη.


The eigenvectors of the Hamiltonian are derived by diagonalizing the transfer matrix, which is achieved by the Bethe-ansatz method <cit.>. The eigenenergies are written in terms of the eigenvalues of the transfer matrix τ(λ),

    E({λ_j}) = (2 sinhξ_+ sinhξ_- coshη  (sinhη)^2N-1)^-1·τ'(0) - (sinhη)^2 + N (coshη)^2/coshη, 
       τ(λ) = (sinh(λ + η))^2Nsinh(2λ + 2η)/sinh(2λ + η)sinh(λ + ξ_+) sinh(λ + ξ_-) ∏_i=1^n sinh(λ - λ_i - η)/sinh(λ + λ_i + η)
       
        + (sinhλ)^2Nsinh(2λ)/sinh(2λ + 
    η)sinh(λ + η - ξ_+) sinh(λ + η - ξ_-) ∏_i=1^n sinh(λ + λ_i + 2η)/sinh(λ - λ_i),

where λ_i solves a set of the Bethe equations,

    ( sinh(λ_j + η)/sinh(λ_j))^2N   = sinh(λ_j - ξ_+ + η) sinh(λ_j - ξ_- + η)/sinh(λ_j + ξ_+) sinh(λ_j + ξ_-)
         ·∏_0ptk=1k ≠ j^n 
        sinh(λ_j - λ_k + η)/sinh(λ_j - λ_k - η)sinh(λ_j + λ_k + 2η)/sinh(λ_j + λ_k),

for j = 1,2,…,n. 

Since the analytic solutions to the Bethe equations (<ref>) are inaccessible, we instead give numerical results for the energy spectra of the XXZ spin chains (Fig. <ref>). As compared from the spectrum of the effective Hamiltonian (<ref>), the sums of the eigenenergies for H^(+) and H^(-) are indeed embedded in its full spectrum. 

Unlike the rSGA-induced partial solvability discussed in the previous subsection, neither equally-spaced spectrum nor pure-imaginary eigenvalue is observed in the solvable subspace of the XXC model.

This indicates that no oscillating mode exists for the XXC model coupled to the boundary dissipators. The only non-decaying mode is the steady state, corresponding to the zero eigenvalue of the effective Hamiltonian. 

Degenerate steady states are observed for the entire Liouvillian, but the integrable subspace has the unique steady state,

    ρ_ ss^XXC = |11 … 1 ⟩⟨ 11 … 1|,

which is a product state. This is also a completely separated state from the other states by the HSF, and therefore initial states never reach the integrable steady state unless it is the steady state itself. 





0
Our main focus in this subsection is how partial solvability of the system affects the relaxation phenomena of this non-equilibrium system. 

The first difference shows up in behaviors of the eigenmodes of the Liouvilian. As the Liouvillian for the XXC spin chain coupled to the boundary dissipators (<ref>) is not fully solvable, a generic eigenmode exponentially decays by time (see Figure <ref>). 
On the other hand, the eigenmodes in the solvable subspace persistently survive, some of which are the steady states, while the others are the persistently oscillating modes, which are known to be in the context of another partially solvable Liouvillian exhibiting “the quantum synchronization" <cit.>. 
In the solvable subspace of the Liouvillian, i.e. equivalently the effective Hamiltonian (<ref>), the dissipation terms always act as zero, and therefore, the dissipators are irrelevant in the solvable subspace. That is, any state in this subspace is a so-called “dark state" <cit.>, which does not feel the dissipators. 

Now let us explicitly construct the steady states and oscillating modes consisting of the solvable eigenmodes. 
As the model has an exact map to the XXZ model with diagonal boundaries in each subspace specified by a certain irreducible strings, the solvable eigenmodes are written as the Bethe states <cit.>. We write them as |E_j ⟩⊗ |E_k ⟩, corresponding to each eigenenergy iE({λ_j}) - iE({λ_k}) given in (<ref>). Note that eigenenergies for the effective Hamiltonian (<ref>) are purely imaginary in the solvable subspace. 
It is straightforward that any state in the form of |E_j ⟩⊗ |E_j ⟩ in the doubled Hilbert space, i.e. a pure energy eigenvector |E_j ⟩⟨ E_j| in the original Hilbert space, is a steady state 

    H |E_j ⟩⊗ |E_j ⟩ = (iE({λ_j}) - iE({λ_j}))  |E_j ⟩⊗ |E_j ⟩ = 0 
       ⇔ℒ(|E_j ⟩⟨ E_j|) = 0,

since E({λ_j}) is purely imaginary. 
The same observation can be applied to the diagonal density matrix consisting of the Bethe vectors, which is also the steady state. Suppose p_j > 0 for all j such that ∑_j p_j = 1. Then one obtains 

    H∑_j p_j |E_j ⟩⊗ |E_j ⟩ = ∑_j p_j (iE({λ_j}) - iE({λ_j}))  |E_j ⟩⊗ |E_j ⟩ = 0 
       ⇔ℒ(∑_j p_j |E_j ⟩⟨ E_j|) = 0.


The oscillating modes are, on the other hand, given by a off-diagonal density matrix consisting of the Bethe vectors. For instance, a density matrix consisting of two energy eigenvectors

    H( p_jk |E_j ⟩⊗ |E_k ⟩ + p^*_jk |E_k ⟩⊗ |E_j ⟩) 
       = p_jk( iE({λ_j}) - iE({λ_k}) )  |E_j ⟩⊗ |E_k ⟩ 
        + p^*_jk( iE({λ_k}) - iE({λ_j}) )  |E_k ⟩⊗ |E_j ⟩

exhibits persistent oscillation for p_jk∈ℝ. 




0
As the Liouvillian (<ref>) is a non-trace preserving map, a density matrix is expected to vanish after a long time t →∞. For a generic pure initial state prepared as the superposition of energy eigenstates

    **,

the Liouvillian acts as 

    ***

unless the state vanishes by the measurement operators. 
Indeed, our numerical test tells that the system prepared in an arbitrary initial state decays by time (see Figure <ref>). 

We found, however, there exist a family of density matrix which does not vanish by applying the Liouvillian (<ref>) repeatedly. 
One of such examples is the pure energy eigenstates of the XXC model in the alternating irreducible subspace

    ρ_ ss = |E_n ⟩⟨ E_n|.

This density matrix indeed is the root of the Liouvillian (<ref>), since it commutes with the Hamiltonian and vanishes under the measurement terms. 
Subsequently, any diagonal density matrix consisting of the energy eigenstates in the same subspace

    ρ_ss = ∑_n p_n |E_n ⟩⟨ E_n|, 
       ∑_n p_n = 1,     ∀ n,  p_n ≥ 0,

also becomes the root of the Liouvillian. Thus, the diagonal density matrix of the form (<ref>) is the steady state, which never decays by time, although the Liouvillian is not a trace-preserving map. 

On the other hand, if one prepares the initial state as a superposition of the energy eigenstate in the subspace of the alternating irreducible string

    |ψ⟩ = ∑_n c_n |E_n ⟩,

the pure state |ψ⟩⟨ψ| exhibits oscillations which persistently survives. It is easy to check that this pure state vanishes under the measurement operators, while the commutator with the Hamiltonian survives 

    [H, |ψ⟩⟨ψ|] = ∑_m,n c_m c_n^* (E_m - E_n) |E_m ⟩⟨ E_n|.

Thus, the persistent oscillations are observed for the system prepared in the superposition of a few number of non-degenerate energy eigenstates. 





  §.§.§ Solvable steady states induced by MPO symmetry

In the previous subsection, we have seen that partial solvability emerges for the open XXC model coupled to boundary dissipators, where the eigenmodes including the steady state are exactly calculated. 
If one focuses only on steady states, some of them are analytically derived even if they are out of the solvable subspace. 

The derivation of those “out-of-integrable" steady states can be achieved based on the MPO symmetry we reviewed in Sec. <ref>. 




By following the method originally introduced in <cit.>, we write the steady state in the form of a product of an amplitude operator Ω,

    ρ_ ss = ΩΩ^†.

Then, if the amplitude operator Ω is given by the matrix product form,

    Ω = _a⟨ v_ L| L_a,N… L_a,2 L_a,1 |v_ R⟩_a,

in which the operator L_a,n satisfies the local divergence (<ref>) with the local perturbed XXC Hamiltonian, the amplitude matrix (<ref>) produces the steady-state density matrix, as we see in the following. 

Let us consider the time evolution of the density matrix consisting of the amplitudes matrix (<ref>). Time evolution of any density matrix is described by the Lindblad equation (<ref>). Due to the local divergence relation (<ref>), the commutator term produces only two terms coming from the non-commutativity at the boundaries,

    [H_XXC, Ω]    = _a⟨ v_ L| M_a,N… L_a,2 L_a,1 |v_ R⟩_a 
        - _a⟨ v_ L| L_a,N… L_a,2 M_a,1 |v_ R⟩_a
       = 0,

since the solution exists only for M = 0. 
On the other hand, the boundary dissipators non-trivially act only at the edges, which may cancel the non-commuting terms above,



    _a⟨ v_ L| ⊗_b⟨ v_ L| ( -i M_a,N L_b,N^†_ p + γ_ R,+𝒟_ R,+(L_a,N L_b,N^†_ p) + γ_ R,-𝒟_ R,-(L_a,N L_b,N^†_ p) ) = 0  , , 
       ( i M_a,1 L_b,1^†_ p + γ_ L,+𝒟_ L,+(L_a,1 L_b,1^†_ p) + γ_ L,-𝒟_ L,-(L_a,1 L_b,1^†_ p) ) |v_ R⟩_a ⊗ |v_ R⟩_b = 0.


Here we leave the boundary dissipation rates free, which shall be constrained later for solvability of the steady state. 
The solutions to the divergence relation (<ref>) and the boundary cancellation condition (<ref>) then provide the steady state density matrix in the form (<ref>) and (<ref>). 

We found that both the diagonal and off-diagonal MPO symmetries on the bulk result in the same solution,

    y=0, 
       |x|^2 = γ_ R,+/γ_ R,- |u|^2,     
        |x|^2 = γ_ R,+/γ_ R,- |v|^2, 
       γ_ R,+/γ_ R,- = γ_ L,+/γ_ L,- = ω.

The integrable steady state observed in the previous subsection is realized as a special case with u = v = 0 and ω = 1. 






§ CONCLUDING REMARKS



In this paper, we have introduced a new class of partially solvable open quantum systems. The models are new in the sense that their partial solvability is induced by partial solvability of the Hamiltonians. The systems are coupled to dissipators only at the boundaries, unlike other existing partially solvable models <cit.>, in which the dissipators are attached at every site. 
We showed that partial solvability of the system can be robust against the boundary dissipators under two different kinds of conditions. 

Liouvillians in the first type admit the existence of dark states, i.e., the states which do not feel the effect of dissipators, and hence partial solvability of the Hamiltonian survives in the subspace spanned by these dark states. A pure dark state becomes an exactly solvable steady state of the Liouvillians, while an off-diagonal element density matrix consisting of the dark states provides an eigenmode of the Liouvillian. 
As an example, we especially focused on the AKLT-type model coupled to quasiparticle baths at both edges. As a known fact, the AKLT-type model exhibits four degenerate ground states <cit.> due to the boundary spin fractionalization, and a tower of solvable quasiparticle excited states can be constructed on top of each ground state <cit.>. Among these four towers of the degenerate solvable energy eigenstates, we found that a tower of states under a certain choice of the boundary spins is a set of dark states of the Liouvillian. 
One of the remarkable features of this model is persistent oscillations obtained in local observables, when the initial state is prepared to have a large enough overlap with the solvable state of the Liouvillian. This quantum synchronization is brought by the equally-spaced spectrum of the Liouvillian in the solvable subspace consiting of the dark states, which is inherited by the equally-spaced spectrum of the Hamiltonian due to the rSGA structure. A similar phenomenon has been reported in <cit.> for the Liouvillian which also exhibits the rSGA but whose dissipators are coupled to all the sites of the system. 

The second mechanism which makes Liouvillians partially solvable is the HSF. We are especially interested in the case where the Hamiltonian exhibits the HSF which divides the Hilbert space into exponentially-many subspaces, some of which may survive even if the boundary dissipators are introduced. The key property in extending the notion of partial solvability induced by the HSF is the robustness of some subspaces under site-dependent perturbations. 
We have considered the XXC spin chain as an example, which exhibits the HSF due to the presence of IS. We showed that the effect of the boundary quasiparticle baths, which inject and absorb quasiparticles, can be regarded as the partial solvability preserving perturbations through the thermofield double formalism. 
As a result, the effective Hamiltonian in the solvable subspace of the doubled Hilbert space becomes two decouples integrable XXZ spin chains with imaginary boundary magnetic fields corresponding to the boundary dissipators. Thus, any eigenmode of the Liouviilian in this solvable subspace is accessible via the Bethe-ansatz method, which is numerically verified (Fig. <ref>). 

Besides the eigenmodes in the solvable subspace, several more steady states can be exactly derived by using the MPO symmetry of the XXC Hamiltonian, which is characterized by the local divergence relation (<ref>). As the local divergence relation (<ref>) holds for the XXC model only when it is reduced to the frustration-free condition (M=0 in Eq. (<ref>)), the steady states associated with the MPO symmetry are the dark states, which vanish under the action of the dissipators. We found that the steady state in the solvable subspace is included in the class of solvable steady states associated with the MPO symmetry. 

A new class of partially solvable open quantum systems introduced in this paper paves a way to study non-integrable open quantum systems. At the same time, it also proposes several interesting questions. We list possible future works in the following, by focusing on the HSF-induced partially solvable open quantum systems.
The first question is how large one can make an overlap with the solvable subspace of the Liouvillian, when a “physical" initial state, such as the dimer state and Néel state, is prepared.

This would be addressed by following the method introduced in <cit.>, which allows the overlap between the initial state and the energy eigenstates to be expressed by the Tsuchiya determinant, after the thermofield double formalism is applied. 
The second question we are interested in is how the relaxation time differs between the solvable subspace and unsolvable subspaces. Definitely, the Liouvillian restricted in the solvable subspace has a gapless spectrum, as it matches the energy spectrum of the XXZ model, and this may lead to a non-trivial (non-exponential) relaxation behavior. 
The third question is whether the Kardar-Parisi-Zhang (KPZ) universality class is observed for open quantum systems. As the XXC chain coupled with boundary dissipators is mapped to the two XXZ spin chains in the solvable subspace, it is likely that the KPZ universality class observed for the XXZ spin chain <cit.> emerges also for this open quantum system, if it is robust against the boundary conditions. 





§ ACKNOWLEDGEMENTS

C. M. thanks H. Katsura, C. Paletta, and B. Pozsgay for helpful discussions, on which the idea of the HSF part of this work is mainly based. C. M. acknowledges financial support from JSPS KAKENHI Grant Number JP23K03244. 
N. T. acknowledges support from JST FOREST (Grant No. JP-MJFR2131) 
and JSPS KAKENHI (Grant No. JP24H00191).






0
The steady state condition ℒ(ρ_ ss) = 0 then requires that the operator L satisfies the local divergence relation (<ref>) and the boundary cancellation conditions

    _a⟨ v_ L| ( -i M_a,N L_a,N^†_ phys + 𝒟_ R,+(L_a,N L_a,N^†_ phys) ) = 0  , , 
       ( i M_a,1 L_a,1^†_ phys + 𝒟_ R,-(L_a,1 L_a,1^†_ phys) ) |v_ R⟩_a = 0.

The local divergence relation is solved by the diagonal or non-diagonal L-operator given in (<ref>), which also solves the boundary cancellation conditions if the parameters are set as in Table <ref>. 

0










§ PROOF OF EQ. (<REF>)
 
In this Appendix, we show that any state in the subspace W^(↑,↓) satisfies the dark-state condition (<ref>) in the presence of the boundary dissipators. 

As the dissipators (<ref>) non-trivially act only on the first and/or the Nth site, the dark-state condition in the present case is written as 

    _a⟨ v_ L| ⊗_b⟨ v_ L| 𝒟_ L(A⃗⊗A⃗^†_p) = 0, 
       𝒟_ R(A⃗⊗A⃗^†_p) |v_ R⟩_a ⊗ |v_ R⟩_b = 0,

in which A⃗⊗A⃗^†_p is the three-by-three matrix with the matrix-valued elements,

    A⃗⊗A⃗^†_p = 
        [ A_0 ⊗ A_0^* A_0 ⊗ A_1^* A_0 ⊗ A_2^*; A_1 ⊗ A_0^* A_1 ⊗ A_1^* A_1 ⊗ A_2^*; A_2 ⊗ A_0^* A_2 ⊗ A_1^* A_2 ⊗ A_2^* ].

The definitions of the matrices A_0, A_1, and A_2 are given in (<ref>). 

Then we have 

    𝒟_ L(A⃗⊗A⃗^†_p) = 
        [      A_2 ⊗ A_2^*                0 -1/2 A_0 ⊗ A_2^*;                0                0 -1/2 A_1 ⊗ A_2^*; -1/2 A_2 ⊗ A_0^* -1/2 A_2 ⊗ A_1^*      A_2 ⊗ A_2^* ],

which indicates that the dissipation terms always include the element A_2 proportional to σ^-. Thus, the dark-state condition (<ref>) is satisfied by choosing the boundary vectors as 

    _a⟨ v_ L| = _a⟨↑|,   _b⟨ v_ L| = _b⟨↑|, 
       |v_ R⟩_a = |↓⟩_a,   
        |v_ R⟩_b = |↓⟩_b.














abbrv




