A note on symmetry reductions of the Lindblad equation: transport in constrained open spin chains

We study quantum transport properties of an open Heisenberg XXZ spin 1/2 chain driven by a pair of Lindblad jump operators satisfying a global `microcanonical' constraint, i.e. conserving the total magnetization. We will show that this system has an additional discrete symmetry which is particular to the Liouvillean description of the problem. Such symmetry reduces the dynamics even more than what would be expected in the standard Hilbert space formalism and establishes existence of multiple steady states. Interestingly, numerical simulations of the XXZ model suggest that a pair of distinct non-equilibrium steady states becomes indistinguishable in the thermodynamic limit, and exhibit sub-diffusive spin transport in the easy-axis regime of anisotropy Delta>1.


Introduction
Simulating interacting quantum many-body systems in-, and in particular, out-ofequilibrium is at the forefront of the current experimental and theoretical research (see e.g. [1], [2]). While it has been recognized that treating large isolated quantum-many body systems (with generic physical properties) poses a formidable task to any method of theoretical description or simulation, it has soon become apparent -first within the community of quantum optics [3] -that often a macroscopic number of observables can be treated as quantum noise, and that the behavior of a few essential physical observables can be extracted in terms of the formalism of open quantum systems [4,5]. From mathematical-physics point of view most notable is perhaps the discovery of Lindblad, Gorini, Kossakowski and Sudarshan [6,7], that within the Markovian approximation any time-development of an open system's (reduced) density operator ρ(t), acting on a N −dimensional Hilbert space H describing the (slow) degrees of freedom of interest, can be described by the quantum Liouville equation in the forṁ where  (1) is the most general one which is (i) local-in-time, (ii) respects the positivity ρ(t) ≥ 0, and (iii) conserves the trace tr ρ(t) ≡ 1, and forms the so-called quantum dynamical semi-group. However, the approach of open-quantum-system has only quite recently been applied to the non-equilibrium problems of condensed matter physics [8,9,10,11,12,13,14,15], such as the controversial quantum transport problem in one-dimension [16,17,18]. The central object of concern in physics applications is the fixed point of dynamical semigroup, ρ ∞ = lim t→∞ ρ(t), or the null-vector of Liouvillian Lρ ∞ = 0 (2) and which shall be termed non-equilibrium steady state (NESS). This state is relevant in the asymptotic long time limit of the evolution of the system, and does not decay despite the dissipative nature of the Lindblad equation. The interesting case, when the fixed point is a pure state ρ ∞ = |ψ ∞ ψ ∞ | is of particular appeal in quantum optics and quantum information theory [19], where such a state |ψ ∞ is called a dark state. Sufficient condition for existence of a dark state is that |ψ ∞ is an eigenstate of H and is annihilated by all jump operators, L m |ψ ∞ = 0. Dark states are also known as examples of decoherence-free states, and are of importance in quantum computing [20]. We note that the Lindblad equation (1) can in general describe the system coupled to several thermal (or chemical, magnetic) baths at different values of the thermodynamic potentials, hence the resulting asymptotic states ρ ∞ can be considered as intrinsically nonequilibrium.
Of central importance is the understanding and control of uniqueness of NESS, for example for applications in quantum memories [21] which could be realized in terms of possibly degenerate (non-unique) NESSs. Another potentially very iteresting application of non-unique (multiple) steady states of Liouvillean quantum dynamics could be in detecting non-equilibrium quantum phase transitions, or regions of criticality [24]. Importantly, multiple fixed points of dynamical semigroups are also a crucial part of the concept of decoherence-free subspaces [20], one of the key ideas for fighting errors in quantum computation. Precise conditions for the uniqueness of NESS, i.e. independence of ρ ∞ of the initial condition ρ(0), has been established a while ago [22,23]. Essentially, it follows from the Evans theorem [22] that NESS is unique if and only if the set of operators {H, L m , m = 1, 2, . . .} generates the entire multiplicative operator algebra B(H). One the other hand, the situations where NESS is not unique have been systematically much less explored. For example, in the case of existence of a discrete symmetry, of either (i) LiouvillianL as a whole, or (ii) of all members of the set {H, L m , m = 1, 2, . . .} one might expect non-unique NESSs classified by the eigenvalues of the symmetry operation.
In this paper we will analyze how the existence of symmetries can allow us to treat the transport problem in a boundary driven open anisotropic Heisenberg XXZ spin 1/2 chain with a micro-canonical constraint, i.e. enforcing the exact conservation of total spin (magnetization). Quite interestingly, our numerical results suggest that the two NESSs corresponding to the zero magnetization sector characterized by the eigenvalue of the parity-like symmetry become indistinguishable in the thermodynamic limit, and that in the easy-axis regime (of anisotropy ∆ > 1) the spin transport becomes sub-diffusive (insulating in the thermodynamic limit), quite different from the results [25,26] for the un-constrained boundary-driven XXZ chain. This seems consistent with a recently claimed co-existence of diffusive and insulating transport in a gapped XXZ model [27]. In the Appendix we provide a general mathematical prescription how in both cases (i,ii) of the previous paragraph the symmetry can be facilitated to reduce the problem and block-diagonalize the matrix of the Liouvillian superoperator, whereas only in the second case (ii) the non-uniqueness of NESS is guaranteed by the existence of a symmetry operator and distinct Liouvillian fixed points can be labelled by the eigenvalues of the symmetry operation.

Symmetric boundary driven open XXZ spin chains
In order to study a simple system which does not have a unique NESS and to provide some additional insight into the debate on quantum transport in one dimension [16,27], we shall discuss a finite open anisotropic Heisenberg XXZ spin 1/2 chain on n sites, with the Hamiltonian σ x,y,z i designate a set of standard Pauli matrices corresponding to physical sites i = 1, 2, . . . , n and acting on a Hilbert space H = (C 2 ) ⊗n of dimension N = 2 n . We couple the chain to magnetic reservoirs only through the boundary spins i = 1 and i = n.
Two cases can be considered: (i) Firstly, we take a set of four Lindblad operators where two real parameters Γ > 0 and µ ∈ [0, 1] designate, respectively, the strength of coupling to a pair of magnetic reservoirs and strength of (non-equilibrium) driving. σ ± i = (σ x i ± iσ y i )/2 are the spin-flip operators. (ii) Secondly, we consider a pair of Lindblad operators again parametrized by two coupling parameters Γ, µ having the same interpretation as above. We note that the Liouvillian flow (1) equipped with (5) strictly conserves the total magnetization M = n i=1 σ z i , at the expense of non-local coupling between the first and the last site. In fact, we can imagine Lindblad operators here as incoherent quantum jumps which transfer spin excitations between i = 1 and i = n sites while conserving total number of spin-excitations (quasi-particles). Alternatively, such a model can be interpreted as a spin ring where hopping is fully coherent on bonds (1, 2), (2, 3) . . . , (n − 1, n) while it is fully dissipative (and asymmetric for µ = 0) on one bond (n, 1), and is a particular case of quantum exclusion process complementing the one presented in Ref. [28] (see also [29,30]).
We note that the quantum transport models described by the Hamiltonian (3) and (4), or (5), can be considered as out of equilibrium steady state models with minimal incoherent input, namely with incoherent processes taking place only on the boundary, or on the single bond.
Let P be a permutation operator which exchanges the site i with n − i + 1 for all i. In fact, P can be written as a unitary operator over H with explicit representation in the computational basis (eigenbasis of σ z i , σ z i |a 1 , a 2 . . . a n = (1 − 2a i )|a 1 , a 2 . . . a n for a i ∈ {0, 1}) as P = (a 1 ,a 2 ,...,an)∈{0,1} n |a 1 , a 2 , . . . , a n a n , a n−1 , . . . , a 1 | and can be interpreted as a parity operator. Combining it with spin flips on all sites we define another unitary (parity-like) operator S as Clearly, S is a Z 2 symmetry of the Hamiltonian, as [H, S] = 0 and S 2 = 1.
In case (i) with (4) one notes that SL w 1 S † ≡ŜL w 1 = L w 4 andŜL w 2 = L w 3 hence the Liouvillian flow (1) exhibits a weak symmetry in the sense of Appendix, Eq. (11). Precisely the same case of symmetrically boundary driven XXZ chain (3,4) has been discussed in several recent papers [31,32,25,33,34,35] and has been shown to exhibit interesting non-equilibrium trasport properties. In this case, one can use Evans theorem [22] in order to prove that NESS is unique. The number of invariant subspaces is n S = 2, with symmetry S eigenvalues s 1 = +1 and s 2 = −1. Let the complete set of 4 symmetry adapted operator sub-spaces be now suggestively denoted as {s α , s β } ≡ B α,β . As pointed out in the Appendix (also explaining the notation), the symmetryŜ can still be used in order to reduce the Liouvillian to two blocks Let us now focus on case (ii). S now becomes a strong symmetry (Appendix, Eq. (12)) of the flow (1,5) asŜL s m = L s m for m = 1, 2. In this case, however, we have an additional strong (continuous, U (1)) symmetry, generated by magnetization operator We shall say that the flow (1,5) obeys a micro-canonical constraint. Thus we can use the Theorem 4.1 of Appendix in order to prove the existence of a pair of distinct NESSs for each eigenvalue of total magnetization s z ∈ {−n, −n + 2, . . . , n − 2, n}, i.e. we have a 2N + 1 dimensional convex set of fixed points of the Liouvillian flow. ‡ As from the quantum transport point of view the most interesting is the zero magnetization sector, we will in the following fix s z = 0 and assume the size n to be even. There we have a line of degenerate NESS parametrized by a number from a unit interval u ∈ [0, 1]

Numerical results and transport properties
Let us now turn to the physical properties of our microcanonically constrained open XXZ chain, in particular to the steady state spin current and magnetization profiles. The spin current operator corresponding to the bond (i, i + 1) is defined as, and obeys the operator continuity equation of the form d  we used exact numerical diagonalization, while for n = 12, 14 and 16 we used a wavefunction Monte-Carlo approach called the method of quantum trajectories, as outlined in Appendix A of Ref. [32] (see also Ref. [14]). An efficient Trotter expansion of the propagator e iHδt with complex coefficients [36] has been employed. The stochastic timedependent Shrödinger equation has been simulated until the current was equal on all bonds, and the statistical error was small, both within the accuracy better than 1%.
We consider three characteristic values of the anisotropy parameter: ∆ = 1/2, ∆ = 1, and ∆ = 2, where the linear-response-transport in the grand-canonical ensemble (or in an open chain without the microcanonical constraint) has been found to be, respectively, ballistic [34], anomalous (super-diffusive/sub-ballistic) [35], and diffusive [18,26]. On the other hand, linear-response transport in the micro-canonical ensemble in the latter case (∆ = 2) has been suggested to be sub-diffusive (insulating) [37]. We chose the bath-coupling parameters as Γ = 1, µ = 0.2, which put our setup in a nearlinear-response regime, however, we also tried other values of bath parameters and the results did not change qualitatively. By considering a large driving parameter close to the maximum µ ≈ 1, we have found negative differential conductance for both symmetry subspaces (for n = 8), consistent with the behavior found in the unconstrained model In the midle row of plots (d,e,f) we compare, again for ∆ = 1/2, 1, 2, and for the largest system size attainable, the magnetization profiles from the two symmetry sectors (red/blue). [31,32]. The spin chain with four sites is particularly interesting as the subspace {−1, −1} is one-dimensional, i.e, the NESS in this subspace has to be pure state, i.e. a dark state. Being a pure state it cannot support current (J = 0) and has a vanishing magnetization on all 4 sites. Explicitly, the dark state then is of the form and does not depend on any system's parameters, ∆, Γ, µ. We did not find any dark states of our flow (1,5) for n ≥ 6. In Fig. 1 we display the behavior of the steady state current J versus the chain length n for the two distinct steady states ρ (1) ∞ and ρ 2 ∞ in the spaces {+1, +1} and {−1, −1}. We find a general trend of convergence of the two curves which suggest that the current might not depend on the symmetry sector in the thermodynamic limit and would therefore be unique. The same applies to magnetization profiles displayed in Fig. 2. This implies a kind of ergodicity and suggests physical irrelevance of the symmetry S in the thermodynamic limit. Furthermore, in the regime ∆ = 1/2 the current J(n) tends to converge to a constant suggesting a ballistic behavior, for ∆ = 1 J(n) seems to slowly (perhaps slightly super-diffusively) decrease with n, while for ∆ = 2, J(n) decreases fast. As far as our numerics can suggest it seems that in the regime ∆ = 2, the decrease of J(n) is definitely much faster than 1/n which is compatible with the suggested insulating behavior [37]. In Fig. 2 we make a detailed analysis of the magnetization profiles in both symmetry sectors, and for different system sizes. We find, consistently, that the profiles in the regimes ∆ = 1/2, ∆ = 1, and ∆ = 2, display ballistic, (slightly super-)diffusive, and sub-diffusive (insulating) behavior, respectively. In is perhaps interesting to remark that in the isotropic case ∆ = 1, for the largest system size that we could simulate, the profile seems to be almost perfectly linear which might also be compatible with thermodynamically normal diffusive behavior in this regime (perhaps related to [38]).
In Fig. 3 we focus on the dependence of the relaxation rates of ρ α (t) to the respective steady state ρ α ∞ on the chain length n. The relaxation rate is determined as the Liouvillean spectral gap R = min λ j =0 (−Re λ j ) where λ j denote the eigenvalues of the Liouvillian superoperatorL in the respective symmetry subspace ({+1, +1} or {−1, −1}). We find consistently that in the {+1, +1} sector the spectral gap always quickly decays with n, perhaps faster than n −3 , which is the gap-scaling derived analytically or numerically suggested for some integrable spin chains with boundary driving [13,35]. On the other hand, we are unable to conclude anything meaningful on the scaling of the gap in the {−1, −1} symmetry sector, where the behavior of R(n) may not even be monotonic. To complement the spectral information we plot in Fig. 4 the set of all Liouvillian eigenvalues λ j which lie sufficiently close to the imaginary line   (i.e. with sufficiently small damping rates −Reλ j ). Further, we considered the density distribution of Liouvillian spectrum projected onto the real line, in fact a cumulative distribution W (r) = λ j θ(r − Re λ j )/ λ j 1 giving the probability that a randomly picked damping rate is larger than −r. Here θ(r) designates a Heaviside step function. The plot of W (r) for the two symmetry sectors shown in Fig. 5 in the easy-axis regime ∆ = 2 suggests existence of two pseudo-gaps, particularly clearly for the {−1, −1} symmetry sector.

Discussion and conclusion
We have presented a small but potentially very useful observation on the symmetry reduction of the Linbdblad master equations describing Markovian open quantum systems. We classified quantum Liouvillian symmetries as weak or strong, depending, respectively, on whether only the generator of master equation as a whole commutes with the symmetry operation, or whether each of the Lindblad (jump) operators, and the Hamiltonian, commute with the symmetry individually. We have shown that only existence of a strong symmetry implies nontrivial symmetry reductions inside the space of steady states (fixed points of the Markovian semi-group). We have also provided a non-trivial example of our construction in the context of quantum transport problem for the strongly interacting XXZ spin 1/2 chain. Namely defining an open XXZ chain with a micro-canonical constraint, we have shown that it exhibits a strong parity-type symmetry and outlined numerical computation of distinct steady states and the corresponding physical observables. Our numerical results seem to suggest that microcanonically constrained open XXZ chain exhibits a sub-diffusive (thermodynamically insulating) behavior in the gapped Ising-like (easy axis) phase.
We note that both, weak and strong symmetries can provide practical advantage in description of quantum Liouvillian flows, as they both allow to block-diagonalize and hence reduce the dimension of the Liouvillian. One can also have a combination of both, namely independent generators of weak and strong symmetries, or non-abelian symmetries of either type. In particular, the study of non-abelian Liouvillian symmetries should be an important line of future research.

Acknowledgements
TP acknowledges useful discussions with Wojciech de Roeck and Herbert Spohn. The work has been financially supported by the grants P1-0044 and J1-2208 of the Slovenian research agency (ARRS).

Appendix: Symmetry reductions of quantum Liouvillian dynamics
Let us consider two cases of symmetric Lindblad master equations with two different kinds of dissipators of the flow (1): (i) For the first kind, which we will call a weak symmetry, we assume that there exist a unitary operator S over H, such that for any x ∈ B(H).
(ii) For the second kind, which we will call a strong symmetry, we assume that a Clearly (ii) implies (i), but not vice versa. We continue by fixing some notation. Let s α = e iθα , α = 1, 2 . . . n S denote distinct eigenvalues of S, and let H α be the corresponding mutually orthogonal eigenspaces, so that we have a complete symmetry decomposition of the Hilbert space The adjoint representation of S on Hilbert space B(H) shall be denoted asŜ, The spectrum of a super-operatorŜ is formed of all possible products s αsβ = e i(θα−θ β ) , sinceŜ(|ψ φ|) = s αsβ |ψ φ| for any |ψ ∈ H α , |φ ∈ H β . Let s ν , ν = 1, 2, . . . nŜ denote the distinct eigenvalues ofŜ, denoted such that s 1 = 1, and s ν = 1 for all ν = 1. Note that the following bounds always hold The lower bound nŜ = n S is reached when {s α } are roots of unity, i.e. in the case of Z n S symmetry obeying S n S = 1. The upper bound nŜ = n S (n S − 1) is reached when there is no arithmetic structure in the spectrum {s α } which may happen when S generates a continuous group of symmetry transformations. Let be the corresponding symmetry decomposition of the operator space B(H), where thê S-eigenspaces B ν are mutually orthogonal in the Hilbert-Schmidt sense. Clearly, having a weak symmetry (11) is equivalent to having commuting superoperatorsLŜ =ŜL (17) which means that the Liouvillian is block diagonalized with respect to the decomposition (16), namelyL where B α,β = {|ψ φ|; |ψ ∈ H α , |φ ∈ H β } for α, β = 1, 2, . . . n S .
(ii) We have at least n S distinct fixed points (NESSs), which we can label by s α , namely each diagonal operator space contains at least one fixed point, ρ α ∞ ∈ B α,α for α = 1, 2, . . . n S .
To prove (ii) we note again that operators with non-vanishing trace can only be contained in diagonal spaces B α,α , due to mutual orthogonality of H α . Hence starting from some ρ α (0) ∈ B α,α with tr ρ α (0) = 1 (say ρ(0) = P α / tr P α where P α is orthogonal projector to the subspace H α ), the flow (1) yields a fixed point ρ α ∞ in B α,α for any α = 1, . . . , n S . Again according to Evans theorem [22] this fixed point can be unique, or further degenerate. In any case, we have at least n S distinct steady states, which may be labelled by the symmetry eigenvalue s α . QED Clearly, the direct sum of all invariant subspaces give the entire operator space and the quotient invariant spaces (16) are given by partial direct sums We remark that (non-positive) initial conditions ρ(0) supported in the non-diagonal invariant spaces B α,β , α = β, which have vanishing trace, could in principle provide extra degeneracy of the null space ofL, i.e. an extra dimension of the convex set of NESSs, or simply complete the description of the relaxation process for a general initial condition ρ(0). In other words, we do not see any argument which would forbid the spectra of block-operatorsL| B α,β to contain 0 even for α = β, although this should be -if possible at all -an exceptional (non-generic) situation.