Decoherence and relaxation of topological states in extended quantum Ising models

We study the decoherence and the relaxation dynamics of topological states in an extended class of quantum Ising chains which can present a manyfold ground state subspace. The leading interaction of the spins with the environment is assumed to be the local fluctuations of the transverse magnetic field. By deriving the Lindblad equation using the many-body states, we investigate the relation between decoherence, energy relaxation and topology. In particular, in the topological phase and at low temperature, we analyze the dephasing rates between the different degenerate ground states.


Introduction
The quantum Ising chains introduced in quantum magnetism [1][2][3][4][5][6] represent a class of exactly solvable many-body systems [7] that exemplifies one-dimensional quantum phase transitions [8][9][10][11][12]. More recently, the quantum Ising model was studied in the non-equilibrium regime to investigate the dynamical behavior of quantum phase transitions, e.g. the quenching in a driven Ising chain [13][14][15][16][17][18], the Kibble-Zurek mechanism [19,20], the Loschmidt echo of a single impurity coupled to the Ising chain [21], the engineered quantum transfer [22], the quantum superposition of topological defects [23], the decoherence dynamics in the strong coupling regime [24] as well as the role of quantum correlations in quantum phase transitions [25][26][27]. Importantly, the generalized class of Ising models can be characterized by a topological number [28][29][30][31][32] and, in the topologically nontrivial phase, localized states can occur at the end of an open chain [1,4] or at the interface separating regions with different topological number [33]. This is associated to the ground state degeneracy in the limit of long chains. For instance, in the case of the XY Ising chain, these end-states correspond to the Majorana zero mode of the one-dimensional fermionic Kitaev model [34][35][36]. Depending on the topological number, the extended models of Ising chains can present more than two end-states, viz. several Majorana zero modes [28-32, 37, 38]. The topology can even change by simply adding a single impurity at one end of an open chain [39].
The quantum Ising model has been experimentally implemented using neutral atoms in optical lattices [40], lattices of trapped ions [41][42][43], Rydberg atoms [44], Josephson junctions [45] and superconducting qubits [46][47][48]. These realizations have to be considered as open quantum systems [49,50] since the degree of freedom corresponding to the spin can be readily affected by interaction with the environment. In general, the interplay between dissipation and interactions in quantum many-body systems presents a rich phenomenology [51][52][53][54][55][56][57]. A quantum reservoir-engineering can lead to desired quantum states [58][59][60][61]. Dynamical instabilities can occur in the phase diagram of driven dissipative systems [62][63][64][65]. The decoherence and relaxation dynamics can be characterized by a slow, algebraic decay [66] or by anomalous diffusion [67]. Other interesting effects are the formation of maximally entangled states protected against phase-flip noise [68] and the non-monotonic critical line in the phase diagram of a system with competing dissipative interactions [69].
Although, a priori, spin lattices synthesized in mesoscopic devices can encode Majorana states [70][71][72], which have potentially application in topological quantum computation, the dissipative interaction affecting such systems distinguish them from other realizations. In topological insulators and semiconducting nanowires, Majorana states are protected by fermion parity conservation against the dephasing induced by bosonic fluctuations. By contrast, when one transforms the Ising chains into the fermionic lattices via the Jordan-Wigner transformation, one has to transform consistently the spin operators coupled to the environment. As a consequence, the system is not anymore topologically protected and the dissipation in the transformed model can induce, for instance, inelastic transitions between states of different parity.
Remarkably, the parity still plays a crucial role when the spins have a longitudinal dissipative coupling, namely they are coupled to the environment via the same spin component coupled to the transverse magnetic field, see Fig. 1. This model of longitudinal dissipation was considered to address the quantum phase transition using the mean field approximation [73] or an exact approach based on path integral [74]. This kind of dissipation was also analyzed to investigate dynamical phase transitions [75][76][77][78], non-equilibrium states in presence of temperature differences [79], quantum diffusion [80] and relaxation dynamics in the strong coupling limit [24].
For a single spin/qubit, the dissipation in the Markovian regime can be characterized by two time constants: the energy relaxation time T 1 (the characteristic time scale in which the spin releases energy to the environment) and the dephasing time T 2 (the time scale after which a coherent superposition of two quantum states reduces to a statistical mixture). In this work, we analyze the energy relaxation and the decoherence dynamics for an extended class of quantum Ising chains formed by N spins under the effect of longitudinal dissipative interaction. This represents the natural dephasing mechanism for the spins in absence of interaction. Such regime can occur when, for instance, the individual energy relaxation time of the single qubit T 1 , is much larger than the individual dephasing time T 2 , with T 1 and T 2 defined in absence of interactions. In general, one has T 2 < T 1 and, in some cases, one can also approach T 2 T 1 , e.g. for flux and fluxonium qubits [81]. Assuming this regime, we consider only the longitudinal coupling as the dominant interaction with the surrounding environment. In order to show that the topological protection is conserved in presence of this longitudinal dissipative interaction, we derive the appropriate Lindblad equation for the many-body system.
In the limit of low temperature, we investigate the correlations between the topology of the spin chain, characterized by a winding number g, and the decoherence in the multidimensional ground state subspace. For the simplest case of the transverse Ising model (g = 1), we discuss the crossover from the trivial to the topological phase. We distinguish different contributions, of thermal or topological origin, appearing in the dephasing rate for an initial state given by a coherent superposition of the ground state and the zero-energy excitation. In the topological regime, these two states are almost degenerate for N 1 and the decoherence rate is set by the overlap of the square modulus of the wave functions of the state localized at the left and the right end of the chain (one Majorana zero mode). This term decreases exponentially by increasing the chain length N such that the coherent superposition survives for a long time, viz. the transient regime to achieve statistical mixture of the two ground states is long compared to the other time scales of the system.
We generalize the results of the simple transverse Ising model by studying an extended model which includes three body, next nearest neighbor interaction, with g = 2 in the topological phase. In this case the ground state subspace is fourfold degenerate, with two ground states in each parity sector (even and odd) and with two zero modes whose wave functions are  (red lines). The three-body interaction (blue lines) with coupling constants J (2) x and J (2) y involves the x-and y-components of two next nearest neighbor spins at position n − 1 and n + 1, and the z-component of the intermediate spin at position n.
localized at the ends of the chain. Hence, the extended model opens the possibility to implement adiabatic quantum computation in a given parity ground state subspace. Naturally, this is impossible in the simple transverse Ising model, as the ground states have different parity, which is conserved by the Hamiltonian. Within each parity subspace of the extended model and at low temperature, we find a formula for the decoherence rate that is proportional to a generalized overlap of wave functions involving both Majorana zero modes. This generalized overlap factor still decreases exponentially with the length of the chain in the limit N 1. Finally, in the extended model, the lowest energy excitations in each parity subspace can relax towards one of the two possible ground states. By preparing the system in one of these excited states, we study the decay rates and the final probability of occupation of the different ground states in a (long) transient regime. We show that the latter quantity is associated to the behavior of the wave functions of the Majorana zero modes and of the single particle spectrum in the topological region g = 2.
This work is organized as follows. In Sec. 2 we recall the class of exactly solvable extended Ising models and their topological characterization. Using the Jordan-Wigner transformation, we map the spin model to the fermions model which we diagonalize using the generalized Bogoliubov transformations. In Sec. 3 we discuss the interaction with the local baths and we derive a Lindblad equation starting from the Bogoliubov operators. The Lindblad operators are associated to the transitions between different many-body states of the system. In Sec. 4 we show the results for the transverse Ising model. Afterwards, in Sec. 5 we discuss the results for the simplest extended quantum Ising model characterized by a fourfold ground state degeneracy for which we analyze the dephasing dynamics of the two ground states of same parity. We summarize our results in Sec. 6.

Extended quantum Ising models and topology
In this section we introduce an extended class of quantum Ising models. We set the notation = k B = 1. These models describe one-dimensional interacting spin chains with equal spins and nearest neighbor interaction as well as next nearest neighbor interaction as shown in Fig. 2. We consider chains of finite length with N spins [82]. The Hamiltonian of these chains reads with and the boundary term with the spin operators σ α n at site n with α = x, y, z obeying the algebra [σ α n , σ with σ ± n = σ x n ± iσ y n and the fermionic (spinless) operators satisfying the anticommutation relations {c n , c † m } = δ nm and {c n , x±y = J (2) x ± J (2) y and the parity operator the resulting fermionic Hamiltonian reads x−y (c n+2 c n + h.c.) , x = 0 has g = 0 for B > J (1) x (trivial phase) and g = 1 for B < J Using the parity projection operators P ± = (1 ± P ) /2 with P + + P − = 1 we project the Hamiltonian onto the subspaces of even and odd parity For a closed chain of finite length, owing to the discrete spatial translation invariance, one defines the unitary transformation in the momentum space as c k ± = e −iπ/4 / √ N n e iπk ± n/N c n , whereas k + (k − ) takes all odd (even) integers between −N and N for the even (odd) subspace such that we can write the Hamiltonian (for one parity subspace) in the following form where we have introduced the pseudo spin representation with s x k = (s + k + s − k )/2, and the effective magnetic field x−y sin (kπ/N ) + J (2) x−y sin (2kπ/N ) , x+y cos (kπ/N ) + J In the long chain limit N 1, the Eqs. (11,12) describe a closed curve in the plane B x , B y when one varies parametrically the wave vector k in the first Brillouin zone. This curve is uniquely defined by the set of parameter in the Hamiltonian. The number of closed loops around the origin defines the winding (topological) number of the spin system [28][29][30][31][32]. More specifically, the winding number g in the xz-plane is defined as the line integral on the closed curve spanned parametrically by the vector k and determines the number of clockwise rotations around the origin. Examples are given in Fig. 3. The transverse Ising model with J y = 0 is represented by a circle with radius J (1) x and the center shifted by B, see Fig. 3(a). Hence for J (1) x < B the origin is not within the circle, thus this regime is referred to as the trivial regime with g = 0. For J (1) x > B the origin is within the circle leading to g = 1, thus we are in the topological regime. The topological regime with g = 1 results in a twofold degenerate ground state, which is equivalent to a Majorana zero mode in the fermionic picture for the open chain [28][29][30][31][32]. Other examples with g = 2 are reported in Fig. 3(b) and Fig. 3(c). In these cases, the system has a fourfold degenerate ground state in the open chain with two Majorana zero modes [28][29][30][31][32].
More explicitly, the Hamiltonian of the open chain H c can be expressed in the following diagonal form with the fermionic Bogoliubov operators γ i and the eigenenergies E i . The Bogoliubov operators are determined by the unitary transformation which we define as where the coefficients (wave functions) ψ L/R i,n and the eigenenergies E i are determined by solving numerically the Lieb-Schultz-Mattis equations [1] (see appendix A). In some cases analytic solutions are available. For instance, in the case of the transverse Ising model, the coefficients ψ whereas f N is the normalization constant such that n |ψ L/R i,n | 2 = 1. The energies of the single particle spectrum are given by x cos(κ i ), whereas the possible k values are the solutions of the following transcendental equation: tan(κ i (N + 1)) = J (1) . An example of the single particle spectrum E i of the open transverse Ising chain is shown in Fig. 4(a). The lowest excitation -which hereafter we denote as i = 0 -has imaginary solution for κ 0 = iq 0 in the topological regime (B < J (1) x ) leading to localization of ψ L 0,n and ψ R 0,n at the ends of the chain [1] and with energy E 0 that vanishes E 0 ≈ 0 in the limit N 1. This represents the zero energy mode whose wave function, in the limit of long chain B/J (1) x N , is given by [1] and ψ R 0,n = ψ L 0,N +1−n , namely we have a Majorana zero mode localized at the ends of the chain with 1/q 0 ≈ 1/ ln(J (1) x /B) as decay length. The wave functions ψ L 0,n and ψ R 0,n are plotted in Fig. 5(a). Hence the system has a twofold degenerate ground state with the ground state |GS defined in the trivial regime and the zero-energy Bogoliubov excitation The energy gap is given by E gap = E 0 in the trivial regime B > J (1) x which strongly decreases around J (1) x = B (i.e. the critical value of the quantum phase transition in the thermodynamic limit), see Fig. 4(a). By contrast, in the topological regime B < J (1) x the lowest fermionic excitation E 0 vanishes (twofold degenerate ground state) and the gap is As a second model analyzed in this work we consider the model Hamiltonian introducing an additional interacting term in the transverse Ising model by setting J (2) y = 0). In this case, the system can approach a winding number g = 2, see Fig. 3 x . An example of the single particle spectrum for this case is shown in Fig. 4(b) for a chain with N = 40 spins. Here we observe the appearance of a first zero-energy excitation, which we denote i = 0 1 , close to the point x in which E 0 1 is strongly reduced (this corresponds to the first critical point where the gap closes exactly in the thermodynamic limit). Furthermore a second zero-energy excitation, which we denote i = 0 2 , also appears close to the point B = J (2) x − J (1) x in which E 0 2 is also strongly reduced (this corresponds to the second critical point where the gap closes again in the thermodynamic limit). After this point, in the long chain limit, the ground state subspace is almost fourfold degenerate with the states |GS and γ † 0 1 γ † 0 2 |GS = |0 1 , 0 2 in the even parity subspace and γ † 0 1 |GS = |0 1 and γ † 0 2 |GS = |0 2 in the odd parity subspace. Again, |GS is connected to the single ground state of the trivial regime. In the topological regime B < J (2) x − J (1) x with g = 2, the wave functions of the two zero-energy states localized at the ends of the chain (Majorana zero modes) have the following analytic formulas [28] ψ L and ψ R 0 i ,n = ψ L 0 i ,N +1−n (for i = 1, 2), whereas the coefficients c 1 , c 2 , c 3 are set by the conditions that the wave functions are normalized and orthogonal ( n ψ L/R 0 i ,n (ψ L/R 0 j ,n ) * = δ 0 i ,0 j ). The inverse of the decay lengths associated to the pairs of localized modes are given by An example of the wave functions ψ L 0 1 ,n , ψ R 0 1 ,n and ψ L 0 2 ,n , ψ R 0 2 ,n is plotted in Fig. 5(b). Notice that the wave function associated to the second zero-energy states has an oscillatory behavior with a longer decay length, whereas the first mode has a decay similar to the single zero energy mode of the transverse Ising.

Lindblad equation for the interacting chain
Before discussing the spin chain of interacting spins, we first recall the results for a single spin coupled to a thermal bath which can be expressed as with S α b (hermitian) operators of the bath (to simplify the notation, we neglect the y component). The shifted interaction in the z-component, namely the term σ z −1, does not affect the rest of our analysis. If the bath is an ensemble of independent harmonic oscillators, this shift can be formally removed by unitary (polaron) transformation that displaces the equilibrium position of the oscillators. In the limit of weak coupling with the environment, assuming the Born-Markov approximation and the secular approximation, one can derive the Lindblad equation, see for example [50]. This means that one factorizes the density matrix of the whole system formed by the spin and the bath, with the bath at thermal equilibrium ρ ≈ ρ s ⊗ ρ bath . The relevant quantities in this approach are the Fourier transforms of the correlators of the bath operators at thermal equilibrium Moreover, in the Markov approximation, the memory effects are neglected assuming a fast decay of the bath correlators in comparison to the time scales of the system. In the last step the secular approximation is used, where the fast rotating terms are neglected, as they average out on larger time scales. Using this approach for the single dissipative spin one finds the energy relaxation rate of the single spin which reads 1/T 1 = κ x (2B) + κ x (−2B). The fluctuations of the longitudinal component of the bath operator leads to pure dephasing whose rate is given by Hereafter we assume 1/T φ , the dephasing rate of the single spin, as a given, effective parameter. The total dephasing rate is given by 1/T 2 = (1/T φ ) + 1/(2T 1 ).
In the limit of T φ T 1 , the longitudinal σ z -coupling to the bath is the dominant one. Therefore, if we regard the system at time scale smaller than T 1 , one can neglect the transverse interaction of the qubit with the environment. Hence, we consider the pure dissipative longitudinal interaction affecting the individual spins and the total Hamiltonian reads We remark that, even if this kind of coupling to the environment in σ z -direction leads to a pure dephasing in the non-interacting spin chain, this is not true anymore when we consider the interacting case. In the latter case, the appropriate basis in the perturbative scheme between system and environment are the many-body eigenstates. The latter are not eigenstates, in general, of the local spin operator σ z n . In other words, this interaction can also lead to energy relaxation in the case of interacting spin chains. Note, however, that the dissipative interaction in Eq. (25) commutes with the parity operator and can not induce energy relaxation between states of different parity.
Having in mind long spin chains, we focus on the case in which the local baths are uncorrelated and homogeneous such that we can write which is a realistic assumption for a homogeneous spin chain with locally separated spins with average spacing larger than the correlation length of the fluctuations of the environment.
For open chain conditions, we express the local spin operator σ z n in term of the fermionic Bogoliubov operators with Following the standard approach similar to a single spin [50], using the Markov-Born approximation combined with the secular approximation, one can derive the Lindblad equation.
In the final result, the relevant quantity are the ladder (Lindblad) operators, which can be obtained by considering the spectral decomposition of the coupling operator σ z n − 1 to the local bath at site n. In the interaction picture we write these operators as and the Lindblad operators of the system are The operators C n (ω) rotate with frequency ω in the time evolution of the interaction picture such that the secular approximation can be used: one drops the fast rotating terms. This approximation can be only used when the energy differences in the system are smaller than the relaxation rates of the open system. Notice that, since we assume finite N , this approximation is valid as the energy differences E i − E j of the discrete spectrum remains finite (excluding the zero modes). However one has to treat the zero modes carefully since their energy (for finite N ) becomes exponentially small (nearly degenerate). This implies that, in the case of the topological regime, one can treat the zero modes as effectively degenerate with the ground state in the Lindblad equation, i.e. in Eq. 31 the Lindblad operator with ω = 0 has terms in the sum with E j = E i = 0 for i = 0 1 and j = 0 2 . Ultimately, this ensures the validity of the Lindblad equation for extremely long times in which the interaction picture of the zero modes can be assumed as time-independent (exp(−i(E 0 i + E 0 j )t) ≈ 1), see Eq. 30. The secular approximation becomes inaccurate near the critical points, where we have the transition between the finite energy excitation to the zero energy mode. Thus in the following analysis we exclude the regime near the critical points. In the interaction picture (neglecting the Lamb shift terms), the final Lindblad equation takes the canonical form which reads Inserting Eq. (31) into Eq. (32), one finds the explicit form of the Lindblad equation which is reported in the appendix C. Hereafter we assume Ohmic dissipation for the transverse correlation function with η setting the dissipative coupling strength to the environment and the bosonic function n B = 1/(e ω/Θ − 1) with Θ the temperature. Notice that, in the limit of vanishing frequency and fixed temperature, one has formally 1/T φ = 2ηΘ. In a finite size system, the average energy spacing in the single particle spectrum |E i −E j | scales algebraically with the length N . By contrast, the separation between the energy of zero-energy excitations in the topological phase and the ground state |GS scales exponentially with the length. At such vanishing energy differences, the presence of other source of noise beyond the Ohmic one can become important. Therefore, to take into account this effect, we use a phenomenological approach and we set the "zero frequency" damping rate 1/T φ as an independent parameter. The derivation of the Lindblad equation as given by Eq. (32) (see also appendix C) is one of the main result of this work. This can be applied to any spin chain system which can be diagonalized via the Jordan-Wigner transformation. In the next section we discuss some applications for two specific cases: the transverse Ising model with winding number g = 1 in the topological phase and the extended model with J (2) x > 0 and winding number g = 2 in the topological phase.

Results for the transverse Ising model
We recover the transverse Ising model by setting the parameter J (2) y = 0 in the general class of the extended chain Hamiltonians. We focus on the low temperature limit Θ E gap where the gap is defined as E gap = E 0 in the trivial regime and as E gap = min i E i with i = 0 in the topological regime. Then the occupation of excited states is small and we can restrict to the lowest excitations of the spectrum formed by single or double particle excitation, as schematically shown in Fig. 6.
Hereafter we focus on the decoherence rate for an initial state given by the superposition between the even ground state |GS and the (odd) lowest excitation |0 , which is degenerate to the ground state |GS in the topological regime. Solving the Lindblad equation we find: where the first term arises from the fluctuations of the energy levels and corresponds to "pure dephasing". This is proportional to the overlap of the two wave functions of the zero mode and reads Λ (Ising) The second term in Eq. (34) is associated to thermal fluctuations between the state |0 = γ 0 |GS and the single particle excitations |i = γ † i |GS and reads with energy exchange E i − E 0 , see Fig. 6. Similarly, the third term in Eq. (34) is related to thermal fluctuations between the state |GS and the double particle excitations γ † i γ † j |GS = |i, j or the transitions between the state |0 and |i, j, 0 and energy difference E i + E j , see x < B. Energy transitions can occur between the even ground state |GS and the two particle excitations with energy difference E i + E j where E i or E j can be also E 0 . Energy transitions are also possible between the lowest odd parity state |0 and the single particle excitations with smaller energy difference E i − E 0 or between the single particle excitations and many particle excitations with higher energy difference E i + E j . (b) Schematic view of the spectrum in the transverse Ising open chain in the topological regime J (1) x > B with the zero-energy mode E 0 → 0. Energy transitions can occur between the even ground state |GS and the two particle excitations at Energy transitions are also possible between the (almost) degenerate ground state |0 and the single particle excitation at energy E i − E 0 ≈ E i or with higher excitations at energy E i + E j . Fig. 6. The term Γ d reads Notice that |i, j also include the states with i, j = 0, namely |i, 0 , see Fig. 6. The contribution Γ d is exponentially small at temperature Θ E gap . As shown in Fig. 6(a), Γ d connects transitions between the double particle excitations separated from the ground state by twice the gap energy E i + E j ∼ 2E gap in the trivial phase. Thus the thermal energy is not sufficient to excite the lowest states to these higher excited states since Γ d ∝ exp(−2E gap /Θ). By contrast, Γ s can be relevant even at low temperature in the trivial regime since it involves transition with typical energy difference E i − E 0 E gap , as shown in Fig. 6(a). Thus this interaction is not suppressed for Θ E gap and leads to additional decoherence in the trivial regime. In Fig. 7(a) we plot Γ s that rises with increasing N , since the energy difference between the state |0 to the next excitation |i becomes smaller for larger N , see Fig. 6(a). Finally, the wave functions ψ L 0,n and ψ R 0,n does not correspond to localized states in the trivial regime and hence we have a finite overlap factor of order one Λ (Ising) g=0 = f 2 N N n=1 sin 2 (k 0 n) sin 2 (k 0 (N + 1 − n)) .
This represent a finite pure dephasing contribution to the decoherence rate as plotted in Fig. 7(a). In the topological regime, Γ d has similar behavior as in the trivial regime hence is strongly suppressed as Γ d ≈ exp(−E gap /Θ). Contrary to the previous trivial regime, the rate Γ s is also strongly suppressed by the gap Γ s ≈ exp(−E gap /Θ) since it now involves transition between states |i and the state |0 , the latter almost degenerate with the ground state, Fig. 6(b)). Finally, the pure dephasing contribution Λ (Ising) g is directly related to the overlap of the localized states. In the large N limit, it can be approximated as namely it is exponentially small due to localization of the end states ψ L 0,n and ψ R 0,n . In Fig. 7(b) we plot the pure dephasing term Λ (Ising) g=1 and the sum of the two contributions Γ s + Γ d . Notice the different scale compared to the trivial regime. As expected Λ (Ising) g=1 has a strong dependence on N whereas Γ s + Γ d are almost constant as varying N since their behavior is ruled by the presence of the energy gap.
5 Results for the extended model J (2) x > 0 We discuss now the extended model at finite value J (2) x > 0, as reported in Fig. 3(b),Fig. 4(b) and Fig. 5(b). We focus only on the topological regime with winding number g = 2 and in the zero temperature limit, namely E gap Θ. In this case the system has a fourfold degenerate ground state separated by a gap of order E gap ∼ 2(J (2) x − J (1) x − B). By the analysis of the previous findings, we can restrict the dephasing dynamics in the ground state subspace, as the contribution due to the interaction with excitations scales with exp(−E gap /Θ) and can be neglected in the zero temperature limit. This means the Lindblad operator can be expressed aŝ with P 0 the projector onto the ground state subspace.
In the even subspace we set the populations as p g = GS| ρ |GS and p 0 = 0 1 , 0 2 | ρ |0 1 , 0 2 and the off-diagonal coherent factor ρ g,0 = 0 1 , 0 2 | ρ |GS . Then we obtain the equation and for the coherence factor we have Examples of the behavior of the coefficients Λ + p and Λ + g=2 appearing in time scales of the single qubit dephasing time are reported Fig. 8. Similar expressions are given for the odd subspace in the appendix B. The coefficients Λ + p ,Λ + p and Λ + g=2 are related to the overlap of the Majorana zero modes and they reduce exponentially with the length of the chain N , and with the decay length of the Majorana zero modes Eq. (21). In other words, the dephasing rate of the topological states is exponentially suppressed compared to the dephasing rate of an individual spin 1/T φ . Notice, however, that in the limit t → ∞, the steady state solution of Eq. (41) and Eq. (43) is simply p 0 = p g and ρ g,0 = 0.
In the last part we analyze the relaxation dynamics of the excited subspace. To simplify the notation, we discuss the relaxation in the even subspace and vanishing temperature limit. We set the populations p j,0 1 as the occupation of the excited states in which the excitation j and the zero energy mode 0 1 are occupied. Notice that such states have energy E j in the topological regime. The populations {p j,0 1 } satisfy the following set of coupled rate equations dp j,0 where the rates are given by and the overlapping factor reads χ ± (j,j ) = n ψ R j,n ψ L j ,n ± ψ R j ,n ψ L j,n 2 .
The Eq. (46) describes the relaxation dynamics of the excited states |j, 0 1 which can decay directly towards one of the two ground states |GS or |0 1 , 0 2 or towards one excited state (1) x = 8B and N = 40 (see Fig. 9). The vertical dashed line corresponds to the energy E gap (see text and Fig.4(b)). of lesser energy E j < E j (see Fig. 9). The last term in Eq. (46) is the positive ingoing flux due to the decay of states at energy E j > E j . Eq. (46) is valid in the time scale in which we neglect internal relaxation in the ground state subspace. This is possible since we have a separation of the time scales: the prefactor κ(E j ) in W (j,0 1 )→g and W (j,0 1 )→0 is ruled by the gap κ(E j ) κ(E gap ) and, at the same time, the overlap factor χ ± (j,0 1 ) involves a delocalized, extended state in the chain with one localized state at the end. The overlap factor χ + (j,j ) involves two delocalized states. In other words, there is no exponential suppression of the rate as in the case of the internal dephasing in the ground state subspace. We solved numerically Eq. (46) to obtain the population p j,0 1 (t) with the initial condition p j,0 1 (0) = δ ij , see Fig.9. To complete the description of the relaxation dynamics, we have to write the equations for the populations of the two ground states p g and p 0 One can check that, at long time t ∼ t ∞ with t ∞ 1/W but still t ∞ T φ /Λ + p , the occupations saturate at values p g (t ∞ ) = p 0 (t ∞ ) (with p g (t ∞ ) = 1 − p 0 (t ∞ )) whereas for times t > T φ /Λ + p the occupations of the ground states approach the values p 0 = p g = 1/2. In Fig. 10 we report p g (t ∞ ) for different initial states (i, 0 1 ) of energy E i . The difference between p g and p 0 strongly depends on the initial state for two reasons: (i) different (internal) decay paths towards lower lying excitations E j < E i , (ii) the different overlap with the two ground states. Naturally, higher excited states have more possible ways to decay towards more lower energy excitations which can become relevant and comparable, a priori, to the direct decay channel towards one of the two ground states. To distinguish between the two different mechanisms of dependence on the initial state, we compare the full expression Eq. (51) with the formula p (dir.) g = W (i,0 1 )→g t∞ 0 dt p i,0 1 (t ) which contains only the direct decay from the initial state toward the ground state, see Fig. 10. We observe that the qualitative behavior of p g (t ∞ ) is well reproduced by the p (dir.) g with larger deviations as increasing the energy of the initial state. In particular, the direct decay description captures the different qualitative behavior at low and high energy.
In Fig. 10 p g (t ∞ ) shows a regular behavior as a function of the initial state up to a maximum initial energy E gap after which there are oscillations as increasing the energy of the initial state. This behavior can be explained by observing the single particle spectrum reported in Fig. 4(b): in the topological regime g = 2 (B < J (2) x − J (1) x ) the spectrum has a regular energy spacing up to some energy E gap , reported as dashed line in Fig. 4(b). Above this energy a intersection of two different bundles of excitations appears. As a result of this intersection the wave functions of the excited states have alternating large and small overlaps with the first ground state |GS , leading to the oscillatory behavior in Fig.10 for E i > E gap . This energy E gap plays the role, roughly speaking, of an effective, secondary gap of the system, which closes at the transition for g = 1 → g = 0 (whereas the primary gap with energy E gap closes at the transition g = 2 → g = 1). Hence the topology of the model is not only a ground state property but can also appear in the relaxation dynamics involving (low energy) excited states.

Summary
In order to understand the robustness of the topological properties of spin chains affected by dissipative interaction with the environment, we studied an extended quantum Ising model in which each single spin is subject to a longitudinal dissipative interaction with a local bath.
For the ground state subspace, we derive the formula for the dephasing rates, in a given parity subspace, that incorporate the two Majorana zero modes showing the robustness of the ground state subspace against the transverse dissipative coupling. The inclusion of a more general dissipative coupling (i.e. with longitudinal) is a interesting future perspective. Furthermore we have shown that the topology can also influence the relaxation dynamics of excited states. Here the secondary gap (which appears for g = 2) determines the relaxation behavior of the excited states and the resulting occupation imbalance of the ground states. It would be interesting if similar behavior can be observed in other topological systems.
Although we focus on the extended Ising chain with a specific form of the next nearest neighbor interaction, with winding number g = 2 in the topological phase, our results can be readily extrapolated to understand the relaxation and dephasing dynamics of the whole extended class of quantum Ising chains.

A Diagonalization of the Ising chains
Rewriting the fermionic Hamiltonian in Eq. (6) in matrix representation with t nm = δ n,m 2B − δ n+1,m J (1) and Inserting the transformation of Eq. (15) into the Hamiltonian Eq. (53), we impose that such transformation diagonalizes the Hamiltonian in the form of Eq. (14) and we get the equations and Setting the matrix (T ) nm = t nm and (∆) nm = ∆ nm , and the vectors ( ψ i L,R ) n = ψ L,R i,n the two Eqs. (56,57) are represented in the following matrix form In general case, we have solved numerically the last equations to find the eigenvectors and the respective eigenvalues E i (single particle energy spectrum).
C Explicit formula of the Lindblad equation