Transition between dissipatively stabilized helical states

We analyze a $XXZ$ spin-1/2 chain which is driven dissipatively at its boundaries. The dissipative driving is modelled by Lindblad jump operators which only act on both boundary spins. In the limit of large dissipation, we find that the boundary spins are pinned to a certain value and at special values of the interaction anisotropy, the steady states are formed by a rank-2 mixture of helical states with opposite winding numbers. Contrarily to previous stabilization of topological states, these helical states are not protected by a gap in the spectrum of the Lindbladian. By changing the anisotropy, the transition between these steady states takes place via mixed states of higher rank. In particular, crossing the value of zero anisotropy a totally mixed state is found as the steady state. The transition between the different winding numbers via mixed states can be seen in the light of the transitions between different topological states in dissipatively driven systems. The results are obtained developing a perturbation theory in the inverse dissipative coupling strength and using the numerical exact diagonalization and matrix product state methods.

Over decades, dissipation has been considered as a destructive influence which destroys the coherence properties of quantum systems. Recently, this point of view has been revised, since tailored environments have been employed in order to dissipatively drive a quantum many-body system into a desired steady state, the so-called attractor state [1]. Even if an external perturbation is applied over a certain time window, the system flows back to the attractor state afterwards. Examples of many-body states that can be reached via an attractor dynamics of a tailored environment are Bose-Einstein condensates [2], number squeezed states [3], Tonks-like states [4], superconducting states [5,6], and, more recently, topologically interesting states [7][8][9][10][11]. These comprise Chern insulators [12] and the Hofstadter model of atoms in an optical cavity [9].
Topological states are characterized by the existence of invariants which can only change in steps by a global action on the system. A paradigmatic example is the use of the stepwise change of the electrical resistance in the quantum Hall effect in topological insulators which is employed for the definition of the standard for the electrical resistance [13]. The classification of topological properties in noninteracting closed systems has attracted considerable attention [14][15][16]. In contrast, topological properties in interacting or In open noninteracting systems [8] two important ingredients were identified for reaching stable topologically nontrivial states. The first one is the existence of a dissipative gap, i.e., a gap in the spectrum of the Lindbladian above the steady state. The second one can only be introduced in noninteracting systems and is the so-called purity gap. This gap measures the purity of the most strongly mixed mode of the bulk.
Here, we go far beyond current studies and show how the intriguing interplay of interactions and a tailored dissipative coupling can give access to topologically interesting properties. To do this, we study by exact analytical and numerical methods the paradigmatic spin-1/2 X X Z-quantum spin chain with dissipative boundaries. Previous work has uncovered far-from-equilibrium steady states of a helical nature with remarkable transport properties [17][18][19][20]. In this Rapid Communication we focus on a specific configuration of this system for which the jump operators at the boundary sites are identical and thus lead to an additional reflection symmetry. We find using a perturbative expansion in the limit of large dissipation that rank-2 steady states-formed by helical states with opposite winding numbers-are dissipatively generated at certain discrete values of the anisotropy parameter due to the space reflection symmetry. The winding numbers have integer values and therefore, similar to topological invariants, can only change their values in integer steps. The helical steady states are not protected by a finite gap which is in contrast to topological states in open systems found previously [8]. However, the dissipative attractor dynamics stabilizes these helical states. As one varies the interaction strength a transition between two helical states occurs, which takes place via higher rank mixtures of states to which several different winding numbers contribute. When the anisotropy changes sign, the steady state transits even via a completely mixed state. Our findings rely on both analytical perturbative methods which are valid for any system size and numerical methods for systems up to N = 12 sites.
We describe the X X Z chain with a density operator ρ by the Lindblad master equation Below we shall seth = 1. The first term on the right-hand side describes the unitary evolution due to the X X Z Hamiltonian, Here, S α j = σ α j /2 are the spin-1/2 operators and σ α j the Pauli matrices acting on site j. The parameter is the anisotropy which determines the quantum phases that appear in an isolated system. The identity I is added for convenience. For | | 1 the ground state of the X X Z Hamiltonian is a gapless Tomonaga-Luttinger liquid. For values | | > 1 a gapped phase occurs which corresponds to a ferromagnetic or antiferromagnetic ground state, respectively. N is the number of sites and we assume in the following for convenience N to be an even number.
The second term describes the dissipative coupling to the environment in Lindblad form Here, is the effective dissipation strength, and L j are the jump operators which act only at the boundary sites j = 1 and j = N and target the density matrix belonging to the eigenstate |↑ x of the spin operator in the x direction defined by σ x |↑ x = |↑ x . Explicitly, L 1 = S y 1 + iS z 1 and L N = S y N + iS z N . We can show that in this situation a unique steady state exists [20].
In the Zeno limit of large dissipative coupling → ∞, the boundary spins to lowest order are pinned in the steady state to the states defined by D 1,N [|↑ x 1,N ↑ x | 1,N ] = 0. The dissipation free subspace of the system is thus the whole Hilbert space spanned by the bulk spins and fixed boundary spins 1, N which are collinear and oriented in the positive x direction, i.e., |↑ x .
Previous studies [17,19,21] have found that for many choices of the boundary dissipation a fine tuning of the anisotropy * m = cos[ϕ m + δϕ/(N − 1)] with the angle ϕ m = (2π m)/(N − 1), with m = −N/2 . . . N/2, generates a pure steady state which is a spin-helix state, where δϕ is a twist angle between the targeted boundary polarizations. We use the * to denote the fine-tuned values and the subindex m in order to distinguish the helicity of the arising state. Here, the state on each site is represented in the basis chosen along the z direction and the spin precesses in the XY plane around the z axis. However, this steady state will become unstable if the spin states targeted at the boundaries become collinear δϕ = 0 and m = 0, as is the case in the chosen situation. We found similar results for δϕ = π .
For the situation, where the spins are locked to the dissipation-free subspace, the system can be viewed as a spin chain on a ring, where the site 1 and N are glued to the same site. Within this configuration, important quantities are the winding numbers of the spin along the ring. They can be determined by the discrete Fourier transform where m = −(N − 2)/2, . . . , (N − 2)/2 denotes the winding number around the z axis and the amplitudes w m can be interpreted as the corresponding weights. Due to the symmetry of the considered system, the relation w m = w −m holds. We note that in a finite system states corresponding to different winding numbers can overlap. However, this overlap vanishes exponentially with system size. In the limit of infinite system size the states corresponding to different winding numbers become orthogonal and the winding number corresponds to a topological invariant.
An intriguing behavior can be seen in the von Neumann entropy S = − i p i log 2 (p i ), where p i are the eigenvalues of the density matrix ρ. In Fig. 1 we show the dependence of the von Neumann entropy on the anisotropy for a small system (N = 6) and a strong amplitude of the dissipative driving = 250J. For this small system, exact diagonalization is used to solve the quantum master equation (1). A drastic behavior in the von Neumann entropy can be seen at the values * m = cos ϕ m with the angle ϕ m = 2π m/5 with m ∈ {1, 2}. At * m two amplitudes w ±m of the winding numbers become dominant, whereas the other values become negligible. This signals that a helical state of rank 2 with two opposite winding numbers arises as the steady state.
Below, after showing further numerical analyses, we will use perturbative arguments to analytically derive that the steady state in the Zeno limit is of the form with |s , |a being orthogonal linear combinations of the spin-helix states | ± m | bulk , restricted to sites 2, . . . , N − 1, with opposite chiralities, |s = A s (|m + |−m )| bulk and |a = A a (|m − |−m )| bulk with b, A s , and A a weights. Additional particularities occur at = 1, where the entropy drops to zero, signaling a pure state which is a helical state corresponding to the winding number m = 0, and at = 0 where a totally mixed state appears. This result demonstrates that steady states with different winding numbers can be reached by a fine tuning of the anisotropy. For finite dissipation strength and fine-tuned anisotropies, we find numerically (not shown) that the steady state is approached as tr[ρ 2 ( ) − (ρ (0) ) 2 ] ∝ (J/ ) 2 , where ρ( ) denotes the steady state at a finite value of . The states at * seem not to be protected by a gap in the spectrum of the Lindbladian-defined as the absolute value of the smallest nonvanishing real part of the eigenvalues of the Lindbladianas can be seen from Fig. 2 where we show that the gap in the Zeno limit closes as 1/ . This is in contrast to previous findings, where the topologically interesting states were protected by a gap [8]. However, the stability of the rank-2 helical state follows from the dissipative nature of the Lindbladian dynamics, the so-called attractor dynamics. Since the steady state is unique for the chosen dissipators, any initial quantum state is guaranteed to approach the targeted rank-2 state (6) asymptotically in time. In particular, also if a perturbation is applied over a certain time, the state relaxes back to the rank-2 helical state. Further, the transition from one helical state to the other goes via the intermediate values of the anisotropy. In Fig. 1, this transition is performed via states which are composed of many different winding numbers and have a much larger von Neumann entropy. For the point which is close to * ±2 , we have a relatively slow dependence, whereas the point corresponding to * ±1 has a very steep dependence. Let us note that the behavior around the special points * steepens with increasing system length.
In order to verify that this is not just a particularity of the small system size, we used a purification implementation of the matrix product state (MPS) method for open quantum systems [22][23][24] as described in Ref. [25] to determine the steady states for larger systems. We have chosen to double the system size to N = 12.
To obtain the steady state, we use the time-dependent MPS method based on a second-order Suzuki-Trotter decomposition with time step t to compute the long-time evolution of an arbitrary state which in this case is chosen to be the Néel state. To overcome the problem of slow relaxation during the attractor dynamics we employ a gradual time evolution procedure. As we are only interested in the steady state and the exact dynamics is irrelevant, we first apply an evolution in a fast-relaxing parameter regime to prepare the initial state ρ ini for the final evolution (see Supplemental Material [26] for details).
This enables us to provide simulation results for different parameter ranges of the interaction anisotropy and the dissipative coupling . The simulation is based on an efficient compression scheme that is well controlled by observing the so-called truncation weight. We verified convergence in this parameter and confirmed that our main findings are not affected by the compression. The final time evolution was computed for a duration of T = 1000/J using a maximal truncation error of ε = 10 −12 and a time step t = 0.1/J. The steady-state expectation values of the required observables are extracted by calculating the average over the last 2000 time steps and are shown in Fig. 3.
Also for these larger systems one can nicely see a similar behavior as described for N = 6. As can be seen in Fig. 3(a), the behavior around the point * ±5 shows that only the winding numbers m = ±5 have an appreciable amplitude and the amplitudes of the other winding numbers rise slowly in its neighborhood. This is compatible with the analytically expected rank-2 steady state decomposed of the two different winding numbers. The steepness of the rise of the amplitudes of additional winding numbers at the special points depends on the system size. In particular, with increasing system size the required value of in order to resolve the special point rises. This is accompanied by an exponential increase of the timescales, such that it becomes very difficult to resolve the steady state in the Zeno limit at * for very large system sizes.
The approach of the Zeno limit can be clearly seen in the dependence on the value of . One finds that the expectation value of the boundary spins collapses already for relatively low values of and becomes locked to the expected value of the dissipation-free subspace around the value of ≈ 100J (not shown). This validates the interpretation that in the large limit the system is close to a ring in which the winding numbers can be associated with topological invariants. Further, as shown in Fig. 3(b) for the value * ±5 , the amplitudes of the winding numbers rapidly approach the expected values for the predicted helical state for increasing , i.e., all amplitudes become negligible except for the amplitudes for m ± 5 which remain finite.
In the following we justify analytically the appearance of the steady state of rank 2 occurring in the Zeno limit at the points * . To this end, we expand the density matrix of the steady state in orders of 1/ as ρ( ) = ∞ n=0 ρ (n) −n . Inserting this ansatz into the Lindblad equation, one can decompose the equation in different orders. The zero-order condition leads to the condition that the density matrix of the boundary spins lies in the dissipation free subspace, i.e., In the first order of expansion (see the Appendix), we obtain the condition Here, H eff acts in the Hilbert space of the internal bulk sites 2, . . . , N − 1 only. It is given by a X X Z Hamiltonian with boundary fields For anisotropies * m the helix states | ± m (4), restricted to the internal sites 2 to N − 1, are eigenstates of H eff with eigenvalue 0, i.e., H eff | ± m | bulk = 0. Thus, the condition (7) is fulfilled by the ansatz R 0 = b|s s| + (1 − b)|a a| which has rank 2. We would like to emphasize that the condition (7) derived perturbatively is very useful for locally acting dissipation and might prove to be useful for a variety of different setups.
In order to find the weight b, we investigate the compatibility conditions arising in the second order in 1/ . Among other conditions (see the Appendix) we obtain where the last equality holds for m and N − 1 coprime. The overlap η vanishes exponentially with system size and the predicted rank-2 steady state has contributions of the two helical states | ± m ±m|. Further conditions (see the Appendix) need to be fulfilled by the steady state, so that the rank-2 state Eq. (6) is not necessarily the steady state. Considering our numerical findings (up to N = 13), we come to the conjecture that the state Eq. (6) is the true steady state at the fine-tuned anisotropy * ±m in the Zeno limit, whenever N − 1 and m are coprime. If the coprime condition is violated, i.e., the ratio m/(N − 1) can be simplified, then we find that the nonequilibrium steady state has higher rank r > 2. Details on the occurring states are presented in Ref. [27].
One very interesting open question which remains is what happens to these findings in the thermodynamic limit. In this limit the fine-tuned values of the anisotropy become dense and the states of different winding numbers become close. It would be interesting to see whether the rank-2 steady states remain stable solutions and how a crossing between the different states can take place.
To summarize, we have found that helical states can be the steady states of a X X Z model of finite size which is coupled at its boundaries to dissipation. We see that in this case the helical states are not protected by gaps in the Lindblad spectrum and that the transition between helical states with different winding numbers goes via highly mixed states. This opens the question of whether other examples exist of topologically interesting states in dissipatively driven systems which are not protected by a gap in the Lindbladian. A further interesting direction would be the investigation of these transitions between a topological state in dissipative quantum systems using the quantum Fisher information [28]. In the following we discuss how we can obtain the proposed rank-2 state in Eq. (6) in the main text from these relations. We note that this perturbative argumentation has been used more commonly and typically the derived equations can be solved using the so-called "adiabatic elimination technique" of virtual excitations [30][31][32]. Here, an additional challenge is the large dissipation-free Hilbert space and the complicated structure of the spectrum of the Lindbladian due to the locality of the dissipation. We circumvent this challenge by deriving additional conditions which are simpler to treat. We would like to point out that the derivation here is very general for the cases of local dissipation.
The zeroth-order Eq. (A1) only gives information at the boundary sites and is satisfied by the ansatz . To obtain information about the bulk part of R 0 , we need to consider the higher-order relations. To obtain information from these, it is convenient to decompose the Hamiltonian as an operator acting in the tensor product space H 0 ⊗ H 1 , where H 0 is a Hilbert space of the two boundary spins 1, N, and H 1 is the Hilbert space of the remaining bulk spins 2, . . . , N − 1. We introduce an orthonormal basis e 0 , e 1 , e 2 , e 3 in H 0 by The Hamiltonian with respect to this basis becomes One can show that the matrix elements between the zeroth and third state vanish, i.e., We introduce H 0,0 ≡ H eff which is given by Eq. (8) in the main text. The commutator in Eq. (A2) for k = 0 can be rewritten using this decomposition as Using this representation and taking the trace over the boundary sites the condition simplifies to which is given in Eq. (7) in the main text. The condition can be fulfilled if we assume the form Here, |α are eigenvectors of H eff and ν α are some real valued, non-negative coefficients. They fulfill the condition α ν α = 1 to give Tr[ρ (0) ] = 1. There exist some subtle issues connected to possible degeneracies of H eff . These in particular can lead to the existence of steady states with higher ranks, which goes beyond the scope of the current Rapid Communication [27]. Further, we can use the representation of the commutator in order to obtain information about ρ (1) We obtain where M 1 ⊗ |e 0 e 0 | is an arbitrary element from the kernel of the dissipator D to be determined by higher orders of the recurrence relations. Inserting the above into Eq. (A3) for k = 1, and again using Eq. (A5), we obtain after some algebra Finally, noting H 0,k = H † k,0 (see also Ref. [33] for details), and writing down the matrix elements α|Q|α = 0, we obtain after some straightforward algebra for any value of α, β =α w αβ ν β = ν α β =α w αβ , (A13) w αβ = | β|H 1,0 |α | 2 + | β|H 2,0 |α | 2 .
In Eq. (A13) we recognize the steady-state equation of a Markov process with w αβ being the rate of the transition from the state α to state β. The explicit form of H 1,0 , H 2,0 can be calculated from Eq. (A6) (see, e.g., Ref. [17]) and is given by Note that the index of the spin operators denotes the sites to which the operator is applied. The Perron-Frobenius theorem guarantees an existence of a unique solution of Eq. (A13) with nonnegative entries, which sum up to 1. The quantities ν α thus have the double meaning of the eigenvalues of Eq. (A9) in the original quantum Markov process and of steady-state probabilities of configurations in a classical Markov process with rates w αβ associated with it (see also Ref. [34]). Now, the rank-2 state assumption, in terms of the associated Markov process Eq. (A13), means that the two states 0,1 form a closed set, with weights b, 1 − b which is a generalization of an absorbing state. The closed set property is w 0,β = w 1,β = 0 for all β > 1. We have checked numerically that the closed set property is satisfied for our setup for all N 13, when N − 1 is a prime number [35]. Thus, Eq. (A13) for α = 0, 1 becomes a closed equation for b, i.e., where w αβ = | β|H 1,0 |α | 2 + | β|H 2,0 |α | 2 , from which we obtain the weights b.