Stationary State Degeneracy of Open Quantum Systems with Non-Abelian Symmetries

We study the null space degeneracy of open quantum systems with multiple non-Abelian, strong symmetries. By decomposing the Hilbert space representation of these symmetries into an irreducible representation involving the direct sum of multiple, commuting, invariant subspaces we derive a tight lower bound for the stationary state degeneracy. We apply these results within the context of open quantum many-body systems, presenting three illustrative examples: a fully-connected quantum network, the XXX Heisenberg model and the Hubbard model. We find that the derived bound, which scales at least cubically in the system size the $SU(2)$ symmetric cases, is often saturated. Moreover, our work provides a theory for the systematic block-decomposition of a Liouvillian with non-Abelian symmetries, reducing the computational difficulty involved in diagonalising these objects and exposing a natural, physical structure to the steady states - which we observe in our examples.


Introduction
Understanding ergodicity in quantum many-body systems remains a fundamental task of mathematical physics. The notion of symmetries plays a crucial role in determining the dynamics and long-time behaviour [1,2,3,4,5,6,7]. For closed many-body systems the existence of extensive symmetries (e.g. due to the underlying integrability, or localization of the model) can lead to the absence of both thermalization and ergodicity [8,9,10,11,12,13,14,15].
Understanding non-ergodicity of open quantum systems is seemingly more difficult. For finite-size systems within the Markovian approximation, i.e. when the time dynamics can be described in the Lindblad master equation framework [16,17,18], the question is partially resolved [19,20,21,22]. It is known that specific kinds of symmetries called strong symmetries [20] lead to degeneracy of the (time-independent) stationary state to which the open system evolves into in the long-time limit. This result connecting the structure of the symmetry operator and the degeneracy of the stationary state has allowed research into many interesting questions, e.g. quantum memory storage and manipulation [23], quantum metrology [24], quantum batteries [25], quantum transport [26,27], exact solutions [28], potential probes of molecular structure [29], dissipative state preparation [30,31], correlation functions [32] and thermalization and relaxation [33,34,35]. Furthermore, degeneracy of the stationary states can be a starting point for realizing discrete time crystals under dissipation [36,37,38].
The strong symmetry theorem [20] has thus far been limited to Abelian symmetries. In this article we derive a lower bound for the stationary-state degeneracy of an open quantum system with multiple, non-Abelian strong symmetries. We achieve this through decomposing the Hilbert space representation of these symmetries into a series of irreducible representations. As illustrative examples, we apply our theory to three archetypal models: a quantum network, the open randomly dissipative XXX Heisenberg model and the spin-dephased Fermi-Hubbard model. The models are interesting for understanding transport in molecules and quantum networks [26,39], effects of disorder on relaxation, and superconductivity out-of-equilibrium [40], respectively. We also study the XXX model under collective dissipation as a simple toy example. Using our theory we are able to derive a tight bound for the stationary state degeneracy in each case, finding numerically that this bound is saturated in most situations.
The stationary-state degeneracy in the XXX spin chain and Hubbard model examples, which scales at least cubically in the system size, provides a space for the storage of quantum information -despite the presence of environmental noise. Moreover, our work provides a theory for the systematic block-decomposition of a Liouvillian with non-Abelian symmetries, reducing the computational difficulty involved in diagonalising these objects and exposing a natural, physical structure to the steady states which we observe in our examples.

Theoretical lower bounds on the stationary state degeneracy of an open quantum system
Under the Markov approximation [16,17,41,18], the dynamical evolution of an open quantum system can be described by the Lindblad master equation (here, and in the remainder of this work, we seth = 1) where ρ is the density matrix of the system, L is the Liouvillian superoperator acting on the space of bounded linear operators B(H). The Hamitonian H in Eq. (1) describes the coherent evolution of the system while the Lindblad 'jump' operators L l describe the interactions between the system and the environment with corresponding coupling strengths γ l . We also define the adjoint of the master equation via the standard Heisenberg picture for the observable O O(t) = tr O exp(Lt)ρ(0) = tr exp(L † t)Oρ(0) (expanding the exponential), which also defines the adjoint Liouvillian superatorL † .

Stationary state degeneracy in the presence of Abelian strong symmetries
As a starting point for our theory we describe the nullspace degeneracy of an open quantum system in the presence of a single strong symmetry [20]. A strong symmetry S of the Lindblad equation is a unitary operator satisfying the commutation relations Now suppose S has n s different eigenvalues. Then we can decompose the Hilbert space into the corresponding blocks for each eigenvalue ‡ The Banach space is created by performing a thermo-field doubling and mapping matrix elements |ψ φ| → |ψ ⊗|φ . This 'vectorization' allows for a natural Hilbert-Schmidt inner product to be defined where the Liouvillian is block-diagonal, i.e. L is invariant in each subspace (H α , H β ) which follows from (3). Furthermore, it can be proven, due to trace preservation of the dynamics, that in every diagonal subspace (H α , H α ) there is at least 1 stationary state of L [20]. Hence, the stationary-state degeneracy of a system with a single strong symmetry has a lower bound of n s .

Lower bound on the null space dimension for systems with non-Abelian strong symmetries
We now prove the main result of this article: the lower bound of the null space in the case of multiple, non-Abelian strong symmetries. We start by stating the following simple, but useful, theorem.
where R(S k ) is the representation of the kth strong symmetry in the Banach space §. We then perform a decomposition of the group representation R into a series of irreducible representations, with each representation R i having a dimension of D i (i.e. it is a matrix of size D i ×D i ). The degeneracy (dimension) of the stationary stateLρ ∞ = 0 is then bounded from below by s i=1 D 2 i .
Proof. DefineR i to be an operator on the full Hilbert space acting as R i in the irrep subspace and trivially in the the rest of the Hilbert space. Following the decomposition we have, [H,R i (S k )] = [L l ,R i (S k )] = 0, ∀i, k, l, and by Schur's representation lemma this implies that H, L l , L † l are multiplies of the identity matrix in the R i blocks [42]. This, in turn, implies that where n and m index the elements of the matrixR i (S k ). This means that there are s i=1 D 2 i linearly independent operators that commute with H, L l , L † l . It immediately follows that these operators are in the kernel ofL † , which has the same dimension as the kernel ofL. § In general, due to these symmetries being non-Abelian, we have that [R(S k ), R(S k )] = 0. Remark 1. Note that the individual states in the representation |i, j , with i = 1, . . . , s indexing the irreducible representation and j = 1, . . . , D i indexing the states inside the irreducible representation, are often highly degenerate. By Schur's lemma both H and L l are proportional to the identity matrix on the subspace |i, j . However, in general they can have further symmetry structure within the |i, j subspace (e.g. they can also be multiples of the identity matrix inside these subspaces). In those cases the degeneracy will be larger than the bound; we explore an example of this situation in Sec. 3.2.2. In many situations, however, we find there is no additional structure and that this bound is saturated.
In that case we may dispense with this final requirement of Th. 1.

Remark 3.
Due to the non-Abelian nature of the symmetry group there may exist stationary states ofLρ ∞,k = 0 that are not density matrices, i.e. ρ † ∞,k = ρ ∞,k . However, by nature of the dynamical map in equation (1) every initial density matrix ρ(0) will always remain a valid density matrix and states of the form ρ † ∞,k = ρ ∞,k are components for valid density matrices in the long-time limit. These states are interesting from a quantum information and computing perspective because they represent quantum coherences that can be used to store quantum information [43].
In the following, we illustrate these results with a series of examples.

Examples
Now that we have identified the lower bound in Theorem 1, we apply this result to several qualitative different models: a fully-connected quantum network, the Heisenberg model, and the Hubbard model. Physically, the first model is used to study energy transport in molecules, the second for studying quantum magnetism and the the third one for studying strongly correlated electrons.
In each example, we use exact diagonalization to compute the steady state degeneracy numerically using the QuTiP library [44,45], and compare the numerical results with the lower bound in Theorem 1.

Fully-connected quantum network
As our first example, we consider a fully-connected quantum network used as a simple model for studying the role of symmetries in energy transport in molecules [29,46]. The model describes an exciton hopping between a series of interconnected lattice sites, pictured in figure (1); similar lattice models are frequently used to represent the dynamics of atoms in optical lattices, electrons in quantum-dot arrays, and photons in photonic waveguides or quantum circuits [47]. These types of models have also proven useful in studying excitonic transport in photosynthetic light harvesting systems [48,49]. The Hilbert space is spanned by the basis {|0 , |1 , |2 ,..., |N } where |i denotes the state with an exciton on site i = 1, . . . , N while |0 is the ground state characterised by the absence of any excitation. In this basis, the Hamiltonian of the system reads where ε g , ε and h play the role of ground state, excitation and kinetic energy scales respectively. In order to drive a current through the network some of the lattice sites are connected to thermal baths which each interact with the system via the jump operators These operators describe the absorption and emission of an exciton on site j where µ ± j is the relevant probability of that process. The Liouvillian is invariant under the permutation of any two sites which are not connected to the thermal baths. This symmetry can be described by the exchange operator for sites i and j The exchange operators that don't act on the bath sites (we assume label the bath sites as a, b, ...) then form a series of strong symmetries: Moreover, these operators do not necessarily commute [P i,j , P k,l ] = 0 ∀i, j, k, l and thus constitute a set of non-Abelian strong symmetries in this open quantum system. More generally, every group element A in the permutation group P , which is formed from the different products of exchange elements is a strong symmetry, i.e., These group elements are thus left null eigenvectors of the Liouvillian [20,21] L † (A) = 0.
Now suppose we label the sites which are uncoupled form the baths as 1, 2, ..., n. An arbitrary permutation operator can be block diagonalized as follows where R n (A) is the natural representation of A, and acts on the space spanned by |1 , |2 ,...,|n . Meanwhile, I N −n acts on the remaining space, spanned by |0 , |n + 1 ,|n + 2 ,...,|N . We anticipate this block-diagonal structure will emerge in the stationary states of the model. Moreover, the natural representation is reducible as it can be written as the direct sum of the irreducible, trivial and standard representations whose ranks are 1 and n − 1, respectively. Hence, the permutation operator P is further block diagonalized: We have thus found an irreducible representation of the strong symmetries and since the left kernel of L contains all possible A, by Theorem 1, its dimension must be at least (n − 1) 2 + 1. Furthermore, for every left null vector of L, there is a corresponding right null vector and the dimension of the right kernel is equal to the left kernel dimension. Thus the stationary state degeneracy of this model is at least (n − 1) 2 + 1.
In table 1 we compare this result to those obtained via exact diagonalization of the full Liouvillian. The results are in complete agreement for all the system sizes we are able to reach and the bound is always completely saturated.
Additionally, we can explicitly show the decomposition in equation (18). The nontrivial part of A acts on the space spanned by {|1 , |2 ,...,|n }. To explicitly uncover the decomposition in equation (18) we change from this configuration basis to a new one: where |φ t is the trivial representation and is defined as which satisfies R n (A) |φ t = |φ t . Meanwhile the |ψ i form the standard representation R s (P ) and describe n − 1 linearly independent basis vectors whose coefficients in this same basis sum to 0 and thus φ t |ψ i = 0 ∀i. Note that the exciton is fully delocalised over the sites uncoupled from the baths in all states in this new basis.
In figure 2 we provide example plots of the stationary states of the system following exact diagonalization. The block-diagonal structure of the steady states in the two bases exposes the decompositions in equations (18) and (19), which provide a natural structure to the steady states of the system.
Finally, we note that the Hamiltonian and every Lindblad operators can be likewise block-diagonalized within this decomposition, and take the form: Here, we have set, without loss of generality, the ground state energy ε g to zero.

XXX Heisenberg model
For our second example we turn our attention to a many-body system with continuous non-Abelian symmetries: the quantum Heisenberg model with zero magnetic field. However, most of the following discussions will be valid for general SU (2) symmetric Hamiltonians. We consider the isotropic XXX spin chain in which the spin-spin couplings are the same in all directions where σ x,y,z i denotes the x, y, z Pauli operator on site i. Importantly, the system has an SU (2) symmetry, which can be encovered through the total spin operators where S x,y,z = N i=1 σ x,y,z i is the total spin operator in the corresponding direction. We use 2 to denote the fundamental representation of the SU(2) symmetry on a single site, where, more generally, n denotes a representation of this Lie algebra over an n dimensional space [51]. We then use this notation to decompose the representation of the SU(2) symmetry over an N -site system into a series of irreducible representations. Assuming that the dissipation does not break the SU(2) symmetry of the model, the lower bound of stationary state degeneracy can be determined using Theorem 1. In table 2 we depict this, showing the decomposition of the spin SU(2) symmetry over an N -site system and the corresponding lower bound on the stationary state degeneracy, which grows as O(N 3 ) In order to demonstrate the strength of this bound, we consider coupling the system to two different SU(2) preserving environments.
3.2.1. Random coupling dissipation. In the first case we couple the system to an environment which induces dissipation between arbitrary pairs of sites with random Table 2. Decomposition of the SU(2) symmetry of the N -site XXX Heisenberg model into a series of fundamental, irreducible representations. For example, 2 ⊗2 = 3 ⊕ 1 refers to the decomposition of two spin 1/2s into a spin singlet and triplet. The third column is the corresponding lower bound for the stationary state degeneracy of the open XXX model, provided the dissipation preserves the overall SU(2) symmetry.

N Decomposition of representation
Degeneracy lower bound (complex) strength. Such a model should be interesting for understanding effects of dissipative disorder on thermalization and transport properties in a quantum manybody system. The master equation then takes the forṁ where each M l is an arbitrary complex matrix with i, j indexing its elements. The jump operators L l describe dissipation proportional to the dipole interaction between sites i and j, with M ij l quantifying the strength of this interaction. Whilst not fully realistic, this type of dissipation is instructive to our theory. The Lindblad operator commutes with all generators of the SU(2) symmetry and thus we can identify a series of non-Abelian strong symmetries via the spin operators We test our lower bound by computing the stationary state degeneracy numerically via exact diagonalization. We generate the M ij l by drawing the real and imaginary parts as random numbers with a uniform distribution over the interval [−1, 1]. We then average the stationary state degeneracy over a series of instances of this dissipation. As shown in figure 3, we find that our bound is saturated and provides an exact description of the stationary state degeneracy of the system.

Total spin dissipation.
As a second example of dissipation which preserves the SU(2) symmetry, we consider dissipation induced by the Casimir operator: such that the master equation readṡ Since [S 2 , S x ] = [S 2 , S y ] = [S 2 , S z ] = 0, the SU(2) symmetry is not broken by the dissipation and thus, again, we have a set of non-Abelian strong symmetries. Again, we compute the actual stationary state degeneracy numerically and compare to our theory in table 2. We report the results in figure 3. They indicate that for systems with site number N > 2, the actual number of stationary states is much greater than the lower bound. We can, however, explain this difference within the structure of our theory. As we will discuss now, there is a further structure to the states within each irreducible representation (see Remark 1 in Theorem 1) which, consequently, increases the steady state degeneracy above the lower bound.
Since L = S 2 commutes with the Hamiltonian and L is Hermitian there is a basis in which both H and L are mutually diagonal. Clearly,L is also diagonal in this basis. We can write the Casimir operator as, where j labels states within each irreducible representation and x labels the M i further states in the |i, j subspace, which have well-defined total angular momentum and angular momentum in the z direction. For example, M 1 = 1, M 2 = 4 and M 3 = 5 for the 5-site system whose decomposition is shown in table 2. We note that L is the identity within each |i, j subspace. By Schur's lemma H is proportional to the identity matrix in each subspace i, We can then diagonalize H in the |i, j subspace, where E (i) x are the energies within the irreducible representation i. Note that there may be further accidental degeneracies in the energies, or more generally additional symmetry structure guaranteeing that for some energies E It is clear that |i, j, x is an eigenstate of both H and L with the same respective eigenvalues ∀j and every pair x, x that has the same energy E (i) Hence, assuming that the energies E (i) x are all distinct the stationary state degeneracy is, which is larger than the lower bound of Theorem 1 if one takes into account only the strong SU (2) symmetry and matches the exact-diagonalization results, see figure 3.

Hubbard model
As a final example we consider the Hubbard model under local dissipation. The Hubbard model is a simple but very successful physical model of solid state systems under the tight-binding approximation. For sake of simplicity we focus on the one-dimensional case. The Hamiltonian for an N -site lattice reads where c i,σ and its adjoint are the annihilation and creation operators for a fermion of spin σ ∈ {↑, ↓} on site i. In addition, n i,σ is the number operator for a spin σ on site i. The quantities t and U are, respectively, kinetic and interaction energy scales, with the summation in the kinetic term taken over nearest-neighbour pairs ij on the lattice.
The Hubbard model has a rich symmetry structure which has made it a natural environment for inducing exotic non-equilibrium phases through dissipation [52,53].
The first of these is the spin symmetry, generated by the spin operators which satisfy the SU (2) relations The second, hidden, symmetry is often referred to as the η-symmetry, and relates to spinless quasi-particles (doublons and holons). It plays an important role in the formation of superconducting eigenstates of the Hubbard model [55]. This η-symmetry is generated by the operators which satisfy the relations The two SU(2) symmetries are abelian: To study the representation of this SU(2) × SU(2) symmetry over an N -site lattice, we first need to describe its representation in terms of a single site. Since the two SU(2) symmetries are commuting, we are able to describe their representation with two individual representations. We define a representation of SU(2) × SU(2) as (m, n) where m and n describing the representations of the first and second SU(2) symmetries respectively. For instance, we can decompose the single-site SU(2) × SU(2) group as This means that the 4-dimensional single-site Hilbert space is a direct sum of two invariant subspaces. In the first subspace the representation of the spin symmetry is two dimensional and the representation of η symmetry is the direct sum of two trivial representations, while in the second subspace the converse is true. This can be seen by writing down the spin and η operators on the one-site Hilbert space in explicit matrix form: The basis of the above matrices is {|vac , | ↑↓ , | ↓ , | ↑ }, where the arrows indicate the presence of a fermion of either spin 'up' or 'down'. We now couple the Hubbard Hamitonian in equation (34) to an environment which induces homogeneous, local spin-dephasing. The dynamics is modelled, under the Markov approximation, via the Lindblad master equatioṅ with the Lindblad operators L l = s z l = 1 2 (n l,↑ − n l,↓ ).
A possible experimental implementation for this master equation is discussed in Ref. [53] where the spin dephasing was shown to induce superconducting order in the long-time limit of the system. Under this dephasing, the η symmetry is preserved [L l , η +,−,z ] = 0 whilst the spin SU(2) symmetry is broken into a U(1) symmetry since [L l , S ± ] = 0 and [L l , S z ] = 0. We represent this SU(2) × U(1) symmetry over the one-site space as where the bold numbers denote the representation of the remaining SU(2) symmetry, and the subscript indexes the value of S z (the U(1) symmetry) in that representation. With this notation, we can then construct the representation over the full Hilbert space for N sites. For example, the representation of the 2-site system is which can immediately be extended to N sites. Using this representation we can apply Th. 1 and predict the lower bound of the stationary space dimension. We find that the lower bound of the null space dimension for an N site system is  In table 3, we compare this lower bound with numerical exact-diagonalization results. We find that the bound is exactly saturated for the small systems where calculations are accessible; larger systems are outside numerical tractability. We notice that, for this example, the growth of this bound is O(N 4 ), which is faster than the usual O(N 3 ) growth because of the additional U (1) symmetry.

Conclusion
We have provided a lower bound on the stationary-state degeneracy of open quantum systems with multiple non-Abelian symmetries. This bound is based on the decomposition of the symmetry group into a series of irreducible representations. As examples we have studied an open quantum network, the dephased XXX spin chain and a spin-dephased Hubbard model. By comparing with exact-diagonalization calculations, we have shown that the bound is saturated in most of the examples we looked at, providing crucial evidence of the strength of our bound for open quantum many-body systems. Thus we expect that if additional degeneracies are discovered in the stationary state through e.g. numerical calculations, these may hint at possible hidden symmetries in the model waiting to be discovered. Alternatively, the system may have accidental degeneracies.
Our general results here open many interesting venues for future study. The most pressing is extending the results to non-Abelian dynamical symmetries. Dynamical symmetries in both closed [56], and open quantum many-body systems [52,40,57] have been shown to guarantee absence of any relaxation to stationarity, instead leading to persistent oscillations and complex dynamics. Our work should also be extended beyond the Markovian approximation, to include more widely applicable Redfield master equations [18]. In cases where the open system has only approximate non-Abelian strong symmetries, one may expect metastable stationary states [58,59].
This work also opens up the possibility of applying our symmetry reduction method to the study of dissipative magnetic systems [60], and other two-dimensional systems [6,7], where the high level of symmetry reduction could be beneficial due to a lack of viable computational methods. We also anticipate further work applying our theory to dissipation-induced frustration of open quantum systems [61] due to the high entropy content of the stationary state.
Lastly, our work on the fully-connected quantum network provides a systematic method to study more complex networks with hierarchical structure and topology [46], as well as other molecular structures with a high degree of symmetry [62,29].