Dissipative spin chain as a non-Hermitian Kitaev ladder
Naoyuki Shibata Hosho Katsura








shibata-naoyuki@g.ecc.u-tokyo.ac.jp




 Department of Physics, Graduate School of Science, the University of Tokyo, Hongo, Tokyo 113-0033, Japan






















             


	We derive exact results for the Lindblad equation for a quantum spin chain (one-dimensional quantum compass model) with dephasing noise. The system possesses doubly degenerate nonequilibrium steady states due to the presence of a conserved charge commuting with the Hamiltonian and Lindblad operators. 
We show that the system can be mapped to a non-Hermitian Kitaev model on a two-leg ladder, which is solvable by representing the spins in terms of Majorana fermions. This allows us to study the Liouvillian gap, the inverse of relaxation time, in detail. We find that the Liouvillian gap increases monotonically when the dissipation strength γ is small, while it decreases monotonically for large γ, implying a kind of phase transition in the first decay mode. The Liouvillian gap and the transition point are obtained in closed form in the case where the spin chain is critical. We also obtain the explicit expression for the autocorrelator of the edge spin. The result implies the suppression of decoherence when the spin chain is in the topological regime.











Valid PACS appear here
                             

                              
Dissipative spin chain as a non-Hermitian Kitaev ladder
    Hosho Katsura
    December 30, 2023
=======================================================







§ INTRODUCTION

	With recent advances in quantum engineering, it becomes increasingly important to study how the coupling to the environment affects a system. The time evolution of such an open system can be described by a master equation. 
Under rather general conditions that the evolution is Markovian and completely positive and trace preserving (CPTP), one obtains the Lindblad equation <cit.> for the time-dependent density matrix. In the past, this quantum master equation 
had been mostly used to describe few-particle systems in, e.g., quantum optics. 
However, recent years have witnessed a growing interest in many-particle systems in the Lindblad setting 
<cit.>. 
	
	Although there are several approaches to analyze the Lindblad equation 
such as perturbation theory <cit.> and numerical 
methods <cit.>, 

exact results for the full dynamics is few and far between. 
In some cases <cit.>, the nonequilibrium steady states (NESSs) can be constructed exactly, but it is more challenging to 
completely 
diagonalize the Liouvillian (the generator of the Lindblad equation). 
In this sense, much fewer cases are known as exactly solvable models <cit.>. 
The difficulty 
lies in dealing with the space of linear operators, the dimension of 
which grows more rapidly than that of the Hilbert space. 
This limits the system size 
amenable to exact numerical diagonalization. 
To make matters worse, it is often the case that effective interactions arise from dissipation even when the Hamiltonian itself is reducible to that of a free-particle system. This prevents us from understanding the full dynamics of the system.
	
In this paper, we present an exactly solvable dissipative model which corresponds to the non-Hermitian many-body quantum system. 






The model we propose has a conserved charge that leads to two exact NESSs. Moreover, our model can be seen as a non-Hermitian Kitaev model on a two-leg ladder <cit.>. Therefore, by applying Kitaev's technique <cit.>, our model can be mapped to free Majorana fermions in a static ℤ_2 gauge field, 

which allows us to fully diagonalize the Liouvillian. 

We numerically identify the gauge sectors where the first decay modes live. Then, 
assuming the flux configurations obtained, we derive the exact Liouvillian gap, the inverse of the relaxation time. 
We also study the infinite temperature autocorrelator of an edge spin and obtain its exact formula by applying techniques from combinatorics. 
		


§ MODELS AND NESSS

	We consider the Lindblad equation for the density matrix ρ
	
    ρt=ℒ[ρ]-i[H,ρ]+∑_i(L_iρ L_i^†-12L_i^† L_i,ρ),

where H denotes the Hamiltonian of the one-dimensional quantum compass model <cit.> given by
	
    H=-∑_i=1^N/2 J_xσ_2i-1^xσ_2i^x -∑_i=1^N/2-1 J_yσ_2i^yσ_2i+1^y,

	and L_i=√(γ)σ_i^z  (i=1,…,N) are Lindblad operators. This form of dissipation is known as dephasing noise <cit.> which kills off-diagonal elements of the density matrix, and hence destroys the quantum coherence. Here,

σ^α_j (α=x,y,z) are the Pauli operators at site j, J_x and J_y are the exchange couplings (in the unit of energy), γ≥0 is the dissipation strength parameter, and N is the number of site. We assume 
	N is even and the open boundary conditions are imposed.
It suffices to consider the case J_x,J_y≥ 0, as the other cases can be obtained by an appropriate unitary transformation. 
The operator ℒ[ρ] is called a Liouvillian or a Lindbladian. A NESS is a fixed point of the dynamics Eq. (<ref>), i.e., an eigenstate of the Liouvillian ℒ with eigenvalue 0. 

For our model, there are two steady states ρ_± (1± Q)/2^N, where 1 is an identity matrix and Q∏_iσ_i^z 
	is a conserved charge, i.e.,
	
    [H, Q]=0, [L_i,Q]=[L_i^† ,Q]=0  for ∀ i.


	The proof goes as follows. Because of the Hermiticity of Lindblad operators, there is a trivial NESS, i.e., completely mixed state ρ_c1/2^N. One can easily verify that if ℒ[ρ_c]=0, then ρ = Qρ_c also satisfies ℒ[ρ]=0. 
Although ρ itself is not positive semi-definite which is a necessary condition to be a density operator, one can construct the following operators P_± (1± Q)/2, which are orthogonal projections, and hence positive semi-definite. 
	This then gives (normalized) steady states ρ_± = P_± / P_±. We have checked numerically for small N that they are the unique NESS of the system. 
	



§ MAPPING TO KITAEV LADDER

	A 2^N× 2^N density matrix ρ can be thought of as a 2^2N-dimensional vector (see Appendix <ref> for details). In this sense, we can identify the Liouvillian ℒ for the one-dimensional chain as a non-Hermitian Hamiltonian of a ladder system <cit.> (see Fig. <ref>):
	
    i(ℒ+γ N)≅ H⊗1-1⊗ H+∑_i=1^N iγσ_i^z⊗τ_i^zℋ,

	where the Hilbert space of the RHS is the “Ket⊗Bra space” and τ_i^z is the Pauli matrix for the ith Bra site. The non-unitary terms in the Liouvillian (<ref>) correspond to the non-Hermitian terms in ℋ. 
	
	
		
	We apply to this ladder Hamiltonian (<ref>) the technique by Kitaev <cit.>, which was originally used to solve the quantum spin model on a honeycomb lattice. 
	In order to solve the model, Kitaev developed an elegant technique: substituting Majorana fermion operators for spin operators σ_j^α→ib_j^α c_j. Here, b_j^α and c_j are Majorana operators obeying the Clifford algebra {c_j, c_k}=2δ_jk, {b_j^α, b_k^β}=2δ_jkδ^αβ, and {b_j^α, c_k}=0 with A,B=AB+BA being the anti-commutator. After the mapping, we have a quadratic Hamiltonian of itinerant Majorana fermions 
(c_i's) in each sector specified by the static ℤ_2 gauge field 
	(i.e., each link has the sign ± 1 in the hopping 
	amplitude). 
	Thus, we can diagonalize the Hamiltonian and obtain all eigenvalues and eigenstates, sector by sector. Besides the honeycomb lattice, Kitaev's mapping is applicable to other lattices with a similar Hamiltonian.
	
	Examples include the ladder system <cit.>, which is the Hermitian analog of our model (<ref>).
	
	Next, we define complex fermions f_i,f_i^† each of which is made of two Majorana fermions c_i and d_i: f_i (c_i+id_i)/2,   f_i^† (c_i-id_i)/2. 
	Here, c_i (respectively, d_i) is the Majorana operator for the σ_i (respectively, τ_i) spin. 
Then,
	
	the model is mapped to the Su-Schrieffer-Heeger (SSH) model <cit.> with imaginary chemical potential <cit.>. The Hamiltonian reads (see Appendix <ref> for details)
	
    ℋ(*μ)   =
    		-∑_i=1^Niγμ_i
    		+∑_k,l𝖠_kl f_k^† f_l,

	where *μ=(μ_1,…, μ_N) and 𝖠 is a tridiagonal and complex-symmetric matrix given by
	
    𝖠    2
    		[ iγμ_1   J_x;   J_x iγμ_2   J_y;         J_y iγμ_3     ⋱;           ⋱     ⋱   J_x;         J_x iγμ_N ].

	μ_i's come from the gauge degree of freedom and take the value ± 1.
	The solution of a non-Hermitian quadratic form of fermions is similar to that of a 
	Hermitian one. 
One can construct many-body eigenstates of ℋ(*μ) by just filling single-particle energies of the Hamiltonian, which can be obtained by diagonalizing 𝖠.

Symmetries of the Hamiltonian (<ref>) enable us to restrict the configurations of μ_i's to consider. 
First, ℋ(*μ) and ℋ(-*μ) have the same spectrum because 

the flux configuration in the ladder system is invariant under sending *μ→-*μ. 
(In view of the SSH model, ℋ(*μ) is transformed to ℋ(-*μ) by the charge conjugation f_j→(-1)^j f_j^†.) Second, due to the inversion symmetry, the transformation (μ_1,μ_2,…,μ_N)→(μ_N,μ_N-1,…,μ_1) leaves the spectrum of ℋ unchanged. 
In the following, we only consider the configurations in which the number of positive μ_i's is not less than the number of negative μ_i's. 



§ LIOUVILLIAN GAP

	Let eigenvalues of the Liouvillian ℒ be Λ_i(ℒ). It can be proved <cit.> that all Λ_i(ℒ) satisfy [Λ_i(ℒ)] ≤ 0. A Liouvillian gap g is defined as
	
    g -max_i
    [Λ_i(ℒ)] 0[Λ_i(ℒ)],

	hence, the inverse of the relaxation time. 
	It is clear from Eqs. (<ref>) and (<ref>) that the Liouvillian gap corresponds to the gap between the first and second largest imaginary parts of eigenvalues of ℋ, the former of which is γ N. The configuration which gives the eigenvalue iγ N is μ_i=+1 for all i. 
	
	The reason is as follows. 
	In this configuration, the one-particle energy levels are obtained just by shifting those of the original SSH model by +2iγ. Therefore, we obtain the eigenvalue iγ N by filling all the energy levels. Then, the Liouvillian gap is recast as
	
    g   =γ N-max_*μ
     m(*μ) 0[∑_i=1^N(λ_i+λ_i2-γμ_i)]
       =min_*μ
     m(*μ) 0(2m(*μ)γ-∑_i=1^Nλ_i-λ_i2),

	where λ_i denotes the ith eigenvalue of 𝖠, and m(*μ) the number of μ_i's which are -1.
	
Since at least one of μ_i's must be -1 when m(*μ)  0, we have g ≤ 2γ. 


One might think that we need to consider the case where μ_i=+1 for all i and some single-particle energy levels are empty. However, this is not the case. The Liouvillian gap in this sector must be greater than or equal to g in Eq. (<ref>), as removing one
fermion from an occupied state in this case decreases the imaginary part of the eigenvalue of H by 2γ.
Thus, it suffices to consider configurations different from the one with μ_i=+1 for all i. 

	
	Figure <ref> shows the numerical results of g as a function of γ for various J_y/J_x for a system size N=10. Here, 

“topological,” “critical,” and “trivial” cases refer to the regions J_y/J_x>1, J_y=J_x, and J_y/J_x<1, respectively, in analogy with the 

Hermitian SSH model. 
	From this figure, we can see a kind of phase transition of the first decay mode in every (topological, critical, or trivial) case. We also numerically obtained the chemical-potential configurations which give the first decay mode, as also shown in Fig. <ref>. We do not show all the configurations which give the same eigenvalues. In Appendix <ref>, up to symmetries mentioned above, every configuration which gives the first decay mode is shown.
In the “phase I” of the critical and 
trivial cases, the gap behaves as exactly g=2γ. The reason for this behavior becomes clear in Eq. (<ref>). When γ is small enough, all λ_i satisfy λ_i>0, then g=2γ follows. In the topological case, the situation is 
slightly different. In this case there exists a λ_i which satisfies λ_i<0, and g is smaller than 2γ. For finite N, there is another configuration in the region of γ≪ 1 (see Appendix <ref>), but this region shrinks to zero as N→∞. In the “phase II”, the gap behaves asymptotically g∝ 1/γ in each case. This increase in relaxation time as γ→∞ can be thought of as the Quantum Zeno effect <cit.>.
					

	Our extensive numerical calculation suggests that the chemical-potential configurations which give the first decay modes do not depend on the system size. Then, under this assumption, we can obtain the exact formula for the Liouvillian gap g and the transition point γ_c in the thermodynamic limit of the critical case J_x=J_y=1 (see Appendix <ref> for more details):
	
    g=
    		{[                                                                          2γ                                                                   (0≤γ≤γ_c); 6^1/3(9 γ^2+√(48 γ^6+81γ^4))^2/3-2· 6^2/3γ^23 γ(9 γ^2+√(48 γ^6+81 γ^4))^1/3                                                                     (γ_c≤γ) ]. ,

where
    
    γ_c=√(√(3)-12)≃ 0.605.

We have confirmed that this result agrees well with the numerical one for N=10.



§ AUTOCORRELATOR AT T=∞

	In this section, we study the `infinite temperature' autocorrelator of the edge spin
	
    C_∞(t)⟨σ_1^z(t)σ_1^z(0)⟩_T=∞=12^N(e^tℒ^∗[σ_1^z]σ_1^z),

	where ℒ^∗ is the adjoint operator of the Liouvillian, which describes the time evolution of an operator X as follows:
	
    tX(t)   =ℒ^∗[X(t)]
       i[H,X(t)]+∑_i(L_i^† X(t)L_i-12L_i^† L_i,X(t)).

	In other words, it corresponds to the Heisenberg picture of the open quantum system. A fundamental motivation in quantum engineering is for localized degrees of freedom to maintain coherence over long times. In Ref. <cit.>, this quantity for the transverse-field Ising 

and XYZ chains without dissipation (i.e., closed system) has been studied as the witness of long coherence times for edge spins. The autocorrelator for the dissipative transverse-field Ising model 
was also studied in Ref. <cit.>. However, very few exact results are available for the time evolution of physical quantities <cit.>. Here, we study C_∞(t) for our model with N=∞. Let us consider how the adjoint Liouvillian acts on σ_1^z. For notational simplicity, we set J_x=q/2,  J_y=1/2 and r 4γ. In this case, one finds

    ℒ^∗[P_0]   =q P_1,

and for n≥ 1, 

    ℒ^∗[P_n]   =
    			
    				-qP_n-1-rP_n-P_n+1   (n: odd)
    
    				P_n-1-rP_n+qP_n+1   (n: even)
    			,

where 

    P_0   =σ_1^z,
    
    		P_n   =
    		
    			-σ_1^y(∏_i=2^nσ_i^z)σ_n+1^x    (n: odd)
    
    			-σ_1^y(∏_i=2^nσ_i^z)σ_n+1^y    (n≥ 2  and even)
    		.



















	It is important to note that P_n's are Hermitian and form an orthonormal set, i.e.,
	
    ⟨⟨ P_i,  P_j⟩⟩12^N(P_i^† P_j)=δ_ij.

	The inner product ⟨⟨· ,·⟩⟩ for matrices is called the Hilbert-Schmidt inner product, 
with which Eq. (<ref>) takes the form
	
    C_∞(t)=∑_n=0^∞t^nn!⟨⟨ℒ^∗^n[P_0],  P_0⟩⟩.

	


Now we compute C_∞(t) by considering the so called “Riordan paths” <cit.> (Motzkin paths <cit.> with no horizontal steps at the bottom line) weighted through q and r. (see Fig. <ref>).
To this end, it is useful to consider 

the generating function F(z;q,r) of the weighted Riordan paths, which can be obtained by the so called “Kernel method” <cit.>. 
The autocorrelator in terms of F(z;q,r) reads
	
    C_∞(t)   =∑_n=0^∞t^nn![z^n]F(z;q,r),

	where 
[z^n] f(z) denotes the coefficient of z^n in f(z). 

This can be rewritten as a contour integral 
	
    C_∞(t)   =∑_n=0^∞t^nn!∮ z2πiF(z; q,r)z^n+1
       =12πi∮F(1/w;q,r)we^tw  w,

	where w=1/z. Here we have chosen the contour in the z-plane so that it surrounds the origin and is sufficiently small. As a result, the contour in the w-plane is sufficiently large.  

The final explicit expression for C_∞(t) is cumbersome and is shown in Appendix <ref>. The analytic and finite-size numerical results of C_∞(t) with γ=0.1 in the topological/trivial regime are shown in Fig. <ref>.


From the exact formula, we can derive the inverse of the decay time τ of C_∞(t) (C_∞(t)t→∞∼e^-t/τ) with nonzero r: for 0 < q < 1, one has
	
    τ^-1   =-η_+(q,r),

and for q > 1,
	
    τ^-1   =
    		
    			r   (0<r≤ 1)
    1+r^22r   (1≤ r≤ q+√(q^2-1))
    
    			-η_+(q,r)   (r≥ q+√(q^2-1))
    		,










where η_+(q,r) (-1-r^2+√((1+r^2)^2-4q^2r^2))/(2r). In particular, we find that the decay is suppressed in the topological regime (0<q<1), although C_∞(t) goes to 0 as t→∞ in both regimes. One might note that C_∞ (t) behaves differently from the Liouvillian gap g; g has a cusp for every, i.e., topological, critical, and trivial regime, while C_∞ (t) has non-analytical points only for critical and trivial regime corresponding to q ≥ 1. There is, however, no contradiction between them. This is because g is determined by the NESSs and the first decay mode, while C_∞ (t) is obtained only from the completely mixed state, i.e., the mixture of the two NESSs.



§ CONCLUSIONS

	We have studied the one-dimensional quantum compass model with dephasing noise and obtained the exact steady states. 
We showed that the model can be mapped to the non-Hermitian Kitaev model on a ladder, which is solvable by representing the spins in terms of Majorana fermions. 

This technique allows us to study the Liouvillian gap exactly. In particular, in the critical case where J_x=J_y, the gap in the thermodynamic limit is obtained analytically.    
We have also studied the autocorrelator of the edge spin and obtained 

its exact formula 
using the technique of combinatorics. 


	We would like to thank T. Prosen and M. Žnidarič for valuable comments. H.K. was supported in part by JSPS KAKENHI Grant No. JP18H04478 and No. JP18K03445. N.S. acknowledges support of the Materials Education program for the future leaders in Research, Industry, and Technology (MERIT). 






§ MAPPING OF THE LIOUVILLIAN TO THE NON-HERMITIAN HAMILTONIAN

	In this Appendix, we describe the identification of the Liouvillian with a non-Hermitian Hamiltonian in detail. Let A be a linear operator on Hilbert space ℋ (i.e., A∈End_ℂ(ℋ)). Assuming (ℋ)=𝒩, there is a complete orthonormal basis |ϕ_i⟩_i=1^𝒩 of ℋ, and A can be regarded as an 𝒩×𝒩 matrix
	
    A=∑_i,j=1^𝒩𝖠_ij|ϕ_i⟩⟨ϕ_j|,

	where a sans-serif style 𝖠∈ M_𝒩×𝒩(ℂ) and its element 𝖠_ij is just a c-number. Now, we 

introduce a new space Ket⊗Bra of dimension 𝒩^2 by the following linear map F:
	
    @C=12pt@R=0.2pt
    				F:   End_ℂ(ℋ)[r]   Ket⊗Bra
       90∈   90∈
       |ϕ_i⟩⟨ϕ_j|@|->[r]   |ϕ_i,ϕ_j ⟩⟩

	Note that F depends on the choice of the basis |ϕ_i⟩_i=1^𝒩, but after fixing the basis, F is an isomorphism, i.e., End_ℂ(ℋ) and Ket⊗Bra are in one-to-one correspondence. 

	Let us consider how superoperators (∈End_ℂ(End_ℂ(ℋ))) look like in the Ket⊗Bra space. All superoperators in the Liouvillian ℒ[ρ] have the form
	
    ρ↦ Aρ B.

	This map is rewritten as
	
    ∑_i,j𝖱_ij|ϕ_i⟩⟨ϕ_j|↦∑_i,j(𝖠𝖱𝖡)_ij|ϕ_i⟩⟨ϕ_j|

	for A=∑_i,j𝖠_ij|ϕ_i⟩⟨ϕ_j|, B=∑_i,j𝖡_ij|ϕ_i⟩⟨ϕ_j| and ρ=∑_i,j𝖱_ij|ϕ_i⟩⟨ϕ_j|. In the Ket⊗Bra space, the superoperator can be seen as the following map:
	
		
    ∑_i,j𝖱_ij|ϕ_i,ϕ_j⟩⟩   ↦∑_i,j(𝖠𝖱𝖡)_ij|ϕ_i,ϕ_j⟩⟩
       =[(∑_i,k𝖠_ik|ϕ_i⟩⟩⟨⟨ϕ_k|)_Ket⊗(∑_j,l𝖡_jl^T|ϕ_j⟩⟩⟨⟨ϕ_l|)_Bra]∑_m,n𝖱_mn|ϕ_m,ϕ_n⟩⟩.

	
	Therefore, this superoperator ρ↦ Aρ B can be thought of as the tensor product of two matrix 𝖠⊗𝖡^T. Note that this matrix is basis-dependent because it is not 𝖠⊗𝖡^†. Then, one can identify the Liouvillian ℒ in Eq. (<ref>) as
	
    ℒ   ≅-iH⊗1+i1⊗ H^T
         +∑_i(L_i⊗ L_i^∗-12L_i^† L_i⊗1-121⊗ L_i^TL_i^∗).

	Here, we do not distinguish operators in italic with matrices in sans-serif. For Eq. (<ref>), transpose or conjugate in Eq. (<ref>) does not matter if we choose a basis which diagonalizes σ_i^z's.



§ NON-HERMITIAN QUADRATIC FORM OF FERMIONS

	In a usual closed free-fermion system, we can generally write the Hamiltonian as follows:
	
    H=*f^†𝖠*f,

	where *f=(f_1,…, f_N)^T, *f^†=(f_1^†,…, f_N^†), and 𝖠 is an N× N Hermitian matrix. In general, all eigenvalues of a Hermitian matrix are real, and eigenvectors form an orthonormal basis (after normalization and orthogonalization in degenerate spaces). In other words, 𝖠 can be diagonalized by a unitary matrix 𝖴, and the Hamiltonian is rewritten as
	
    H=*f^†𝖴 diag(a_1,…,a_N) 𝖴^†*f,

	where a_i (i=1,…, N) are the eigenvalues of 𝖠. Then, we can define new operators *f'𝖴^†*f, and it is easily verified that *f' also satisfy anticommutation relations. Therefore, we obtain
	
    H=∑_i=1^N a_i f'_i^† f'_i,

	and all eigenvalues of H are obtained by E=∑_i a_in_i with arbitrary choice of each n_i=0 or  1.
	
	However, some care must be taken when 𝖠 is non-Hermitian. First, eigenvalues of 𝖠 are not necessarily real. Second, eigenvectors with different eigenvalues are in general not orthogonal. Third, left and right eigenvectors with the same eigenvalue are in general not a Hermitian-conjugate to each other. We briefly explain the general prescription for treating non-Hermitian matrices according to Ref. <cit.>. A similar discussion can be found in Ref. <cit.>.

Let 𝖠 be a non-Hermitian and non-degenerate N× N matrix. Remember that if a matrix is non-degenerate, then it is diagonalizable. [We can easily generalize the following discussion to the case where 𝖠 is degenerate but diagonalizable. However, some degenerate matrices cannot be diagonalized, and such situations are sometimes called “coalescence” instead of “degeneracy”. We do not consider “coalescence” here. See, e.g., Ref. <cit.> for more details about coalescence.] Then, 𝖠 has left/right eigenvectors {*l_i^†}_i=1^N and {*r_j}_j=1^N which satisfy
	
    *l_i^†𝖠   =λ_i*l_i^†,
    𝖠*r_j   =ξ_j*r_j.

	(*l_i and *r_j are column vectors.) From these, we obtain
	
    λ_i⟨*l_i, *r_j⟩=*l_i^†𝖠*r_j=ξ_j ⟨*l_i, *r_j⟩
    ∴   (λ_i-ξ_j)⟨*l_i, *r_j⟩=0,

	where ⟨*l_i, *r_j⟩*l_i^†*r_j is the standard inner product of ℂ^N. Because 𝖠 is diagonalizable, each of *l_i_i=1^N and *r_j_j=1^N is a basis of ℂ^N. Therefore, for each i, at least one j satisfies ⟨*l_i, *r_j⟩ 0. Then, we can assume λ_i=ξ_i after relabeling j's. It follows that ⟨*l_i, *r_j⟩=0 if i j and ⟨*l_i, *r_i⟩ 0. Therefore, after “normalization” of *l_i_i=1^N and/or *r_j_j=1^N so as to satisfy ⟨*l_i, *r_i⟩=1, we obtain
	
    ⟨*l_i, *r_j⟩=δ_ij.

	Then, we define
	
    𝖵(*r_1 …*r_N)   (𝖵^-1=
    		[ *l_1^†;      ⋮; *l_N^† ]
    		),

	and 𝖠 is diagonalized as
	
    𝖵^-1𝖠𝖵=diag(λ_1,…, λ_N).

	
	Now, let us return to Eq. (<ref>) with non-Hermitian 𝖠. By diagonalizing 𝖠 as above, we obtain
	
    H=*f^†𝖵 diag(λ_1,…, λ_N) 𝖵^-1*f

	and we define
	
	
    a_i*l_i^†*f,   b_i^†*f^†*r_i,

	then, the Hamiltonian can be written in the form similar to Eq. (<ref>),
	
    H=∑_i=1^N λ_i b_i^† a_i.

	One can easily verify the following anticommutation relations
	
    a_i,b_j^†=δ_ij,  a_i,a_j=b_i^†,b_j^†=0.

	Therefore, a_i and b_i^† can be seen as an annihilation and an creation operator of new fermions, respectively, although a_i^† is not equal to b_i^† unlike Hermitian cases <cit.>. From these anticommutation relations, it follows that b_i^† a_i has eigenvalues 0 and 1, although it is not Hermitian. Then we obtain all eigenvalues of H as E=∑_i λ_in_i with arbitrary choice of each n_i=0 or  1.
	


§ FROM KITAEV LADDER TO SSH MODEL

	By Kitaev's mapping, a non-Hermitian Hamiltonian ℋ can be seen as a model of free Majorana fermions in a static ℤ_2 gauge field. Introducing Majorana fermion operators as σ_j^α→i b_j^α c_j and τ_j^α→ib̃_j^α d_j (i=1,…, N, α=x,y,z), we obtain
	
		
    ℋ   =J_x∑_j=1^N/2[(ib_2j-1^xb_2j^x)(ic_2j-1c_2j)-(ib̃_2j-1^xb̃_2j^x)(id_2j-1d_2j)]+J_y∑_j=1^N/2-1[(ib_2j^yb_2j+1^y)(ic_2jc_2j+1)-(ib̃_2j^yb̃_2j+1^y)(id_2jd_2j+1)]
       -iγ∑_i=1^N(ib_i^zb̃_i^z)(ic_id_i),

	

where the operators, ib_2j-1^xb_2j^x, ib_2j^yb_2j+1^y, ib̃_2j-1^xb̃_2j^x, ib̃_2j^yb̃_2j+1^y, and ib_i^zb̃_i^z commute with the Hamiltonian and their eigenvalues are  ± 1. Therefore the Hilbert space splits into sectors labeled by the eigenvalues of these operators. One can define the flux  through each plaquette by the eigenvalue of the product of b and b̃ operators around it. The spectrum of the Hamiltonian H depends only on the set of fluxes. Thus we can fix all but N-1 signs of the links, as we have N-1 plaquettes in our model. We fix them as
	
    ib_2j-1^xb_2j^x=ib_2j^yb_2j+1^y=-1,  ib̃_2j-1^xb̃_2j^x=ib̃_2j^yb̃_2j+1^y=+1

	and define μ_i=-ib_i^zb̃_i^z to recast ℋ as
	
    ℋ   =-J_x∑_j=1^N/2[(ic_2j-1c_2j)+(id_2j-1d_2j)]
       -J_y∑_j=1^N/2-1[(ic_2jc_2j+1)+(id_2jd_2j+1)]
       +iγ∑_i=1^Nμ_i(ic_id_i).

	Next, we define complex fermions f_i and f_i^†, 

consisting of two Majorana fermions c_i and d_i:
	
    f_ic_i+id_i2,   f_i^†c_i-id_i2.

	It is easy to verify that they satisfy the anticommutation relations
	
    f_i,f_j^†=δ_ij,   f_i, f_j  = f^†_i, f^†_j = 0.

	We then have
	
    ℋ   =2iJ_x∑_i=1^N/2(f_2j^† f_2j-1-f_2j-1^† f_2j)
       +2iJ_y∑_i=1^N/2-1(f_2j+1^† f_2j-f_2j^† f_2j+1)
       +iγ∑_i=1^N μ_i(2f_i^† f_i-1).

	After the unitary transformation f_j→e^i(π/2)jf_j, we obtain Eq. (<ref>).
	


§ THE FIRST DECAY MODES' CONFIGURATIONS

	In Fig. <ref>, we show for each phase only one example of the configurations where the first decay mode lives, but numerical calculation reveals that there are other such configurations. Tab. <ref> shows the numerical results of 

all such configurations up to symmetries mentioned in the main text.
	


§ DERIVATION OF EQ. (<REF>)

	We call the pattern of *μ which satisfies “μ_1=-1, otherwise μ_i=+1” “pattern 1,” and that which satisfies “μ_1=μ_2=-1, otherwise μ_i=+1” “pattern 2.” 𝖠_1 (respectively, 𝖠_2) denotes the matrix 𝖠 with the pattern 1 (respectively, pattern 2). The (unnormalized) eigenstates *v=(v_1,v_2,…, v_N) of 𝖠_1 whose eigenvalues have negative imaginary parts are obtained by the following ansatz
	
    v_2n-1=α^n-1,    v_2n=iβα^n-1   (n=1,…,N/2),

	where α,β∈ℂ and α<1. Letting λ be the eigenvalue of 𝖠_1 corresponding to *v, we obtain the following conditions for this ansatz
	
    left boundary:         -2iγ + 2iJ_xβ=λ_1,
    bulk:         
    			2J_x-2γβ+2J_yα=iβλ_1
    
    			2iJ_yβ+2iγα+2iJ_xαβ=αλ_1
    		.

	Here, we neglect the right boundary condition that is justified in the thermodynamic limit. The solution of λ is
	
    λ_1=-i2γ[-J_y^2+√(8γ^2(2γ^2+J_y^2-2J_x^2)+J_y^4)].

	In the critical case of J_x=J_y=1, λ has negative imaginary part when γ>1/√(2), and the gap for pattern 1 [i.e., the argument in parentheses in Eq. (<ref>) for pattern 1] reads
	
    g_1=
    		
    			2γ   (0≤γ≤ 1/√(2))
    1γ   (1/√(2)≤γ)
    		.

	In a similar way, we obtain the localized solution for pattern 2 by the ansatz
	
    v_1=1, v_2=iδβ, 
    
    			v_2n-1=δα^n-1, v_2n=iδβα^n-1  (n=2,…, N/2)
    		,

	and conditions
	
    left boundary:          
    			-2iγ+2iJ_xδβ=λ_2
    
    		2J_x+2γδβ+2J_yδα=iδβλ_2
    		,
    bulk:          
    			2J_x-2γβ+2J_yα=iβλ_2
    
    			2iJ_yβ+2iγα+2iJ_xαβ=αλ_2
    		.

	From these conditions, we obtain in the critical case the gap for pattern 2 as
	
		
    g_2=
    			
    				4γ   (0≤γ≤√((√(5)-1)/8))
    6^1/3(9 γ^2+√(48 γ^6+81γ^4))^2/3-2· 6^2/3γ^2/3 γ(9 γ^2+√(48 γ^6+81 γ^4))^1/3   (√((√(5)-1)/8)≤γ)
    			.

	
	Finally, we obtain the global gap as
	
		
    g   =min(g_1,g_2)
       =
    			
    				2γ   (0≤γ≤√((√(3)-1)/2))
    6^1/3(9 γ^2+√(48 γ^6+81γ^4))^2/3-2· 6^2/3γ^2/3 γ(9 γ^2+√(48 γ^6+81 γ^4))^1/3   (√((√(3)-1)/2)≤γ)
    			.

	
	Therefore, the transition point γ_c for the critical case is
	
    γ_c=√(√(3)-12)≃ 0.605…

	One may guess that even for topological or trivial case, we can obtain the exact results for the gap in a similar way. However, it would be impossible to obtain the algebraic solutions because in these cases, we need to deal with equations of degree greater than four.
	


§ EXACT FORMULA FOR THE AUTOCORRELATOR WITH DISSIPATION

	The generating function for the weighted Riordan path is obtained by the Kernel method <cit.> as
	
		
    F(z;q,r)-1+(1-q^2+r^2)z^2+√([(1+q^2)z^2+(1+r z)^2]^2-4q^2 z^4)2z[q^2r z^2+(1+r^2)z+r].

	
	Then, the autocorrelator is recast as
	
		
    C_∞(t)   =12πi∮F(1/w; q,r)we^tw  w
       =12πi∮-w^2+1-q^2+r^2+√([(w+r)^2+(1+q)^2][(w+r)^2+(1-q)^2])2r(w-η_+(q,r))(w-η_-(q,r))e^tw  w,

	
	where
	
    η_±(q,r)    -1-r^2±√((1+r^2)^2-4q^2r^2)2r.

	
	The contour of w is chosen as shown in Fig. <ref>, and the final results are
	
		
    C_∞(t)   =
    			-η_+^2+1-q^2+r^2r(η_+-η_-)e^η_+ t
      +e^-r tπ r∫_1-q^1+qf(y,q)(r^-1-r)ycos(yt) +[y^2-(r+η_+)(r+η_-)]sin(yt)[y^2+(r+η_+)^2][y^2+(r+η_-)^2]  y   (0<q≤ 1)
    e^-r tπ r∫_q-1^q+1f(y,q)(r^-1-r)ycos(yt) +[y^2-(r+η_+)(r+η_-)]sin(yt)[y^2+(r+η_+)^2][y^2+(r+η_-)^2]  y   (q≥ 1,  0≤ r<1)
    e^-tcos(√(q^2-1)t)+e^-tπ∫_q-1^q+1f(y,q)y^2-q^2+1sin(yt)  y   (q≥ 1,  r=1)
    -η_+^2+1-q^2+r^2r(η_+-η_-)e^η_+ t--η_-^2+1-q^2+r^2r(η_+-η_-)e^η_- t
      +e^-r tπ r∫_1-q^1+qf(y,q) (r^-1-r)ycos(yt) +[y^2-(r+η_+)(r+η_-)]sin(yt)[y^2+(r+η_+)^2][y^2+(r+η_-)^2]  y   (q≥ 1,  r>1)
    			,

	
	where
	
    f(y,q)=√([(q+1)^2-y^2][y^2-(q-1)^2]).

	In general, it does not have a simpler form. However, in the absence of dissipation, i.e., when r=0, we have
	
		
    C_∞(t;r=0)=
    			
    				1-q^2+1π∫_1-q^1+q√([x^2-(1-q)^2][(1+q)^2-x^2])xcos (xt)   x    (0<q≤ 1)
    1π∫_q-1^q+1√([x^2-(q-1)^2][(q+1)^2-x^2])xcos (xt)   x    (q≥ 1).
    			.

	
	We have confirmed numerically that the integral goes to zero as t→∞ in either case. Therefore, the autocorrelator is non-vanishing as t→∞ in the topological phase, whereas vanishing in the trivial phase. Moreover, when q=1, it takes a simpler form:
	
    C_∞(t;q=1,r=0)=J_1(2t)t,

	where J_1 is the Bessel function of the first kind. From the asymptotic behavior of the Bessel function, we obtain that C_∞ (t) decays as ∼ t^-3/2 for large t.
	


