Stationary state degeneracy of open quantum systems with non-abelian 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 symmetry whose operator representation R(S) satisfies the commutation relations

$$\begin{equation}\left[H,R\left(S\right)\right]=0,\quad \left[{L}_{l},R\left(S\right)\right]=\left[{L}_{l}^{{\dagger}},R\left(S\right)\right]=0,\quad \forall l.\end{equation} \tag{ 3 }$$

Now suppose R(S) has n_s different eigenvalues. Then we can decompose the Hilbert space into the corresponding blocks for each eigenvalue

$$\begin{equation}\mathcal{H}=\underset{\alpha =1}{\overset{{n}_{\mathrm{s}}}{\bigoplus}}{\mathcal{H}}_{\alpha },\end{equation} \tag{ 4 }$$

with α indexing the different blocks. Moreover, we can carry this decomposition through to the space of bounded linear operators $\mathcal{B}\left(\mathcal{H}\right)=\left(\mathcal{H},\mathcal{H}\right)$

$$\begin{equation}\left(\mathcal{H},\mathcal{H}\right)=\underset{\alpha =1}{\overset{{n}_{\mathrm{s}}}{\bigoplus}}\underset{{\alpha }^{\prime }=1}{\overset{{n}_{\mathrm{s}}}{\bigoplus}}\left({\mathcal{H}}_{\alpha },{\mathcal{H}}_{\alpha }^{\prime }\right),\end{equation} \tag{ 5 }$$

where the Liouvillian is block-diagonal, i.e. $\mathcal{L}$ cannot move states between any of the the ${n}_{\mathrm{s}}^{2}$ subspaces which are indexed by α and α'

$$\begin{equation}\mathcal{L}\left({\mathcal{H}}_{\alpha },{\mathcal{H}}_{\alpha }^{\prime }\right)\subseteq \left({\mathcal{H}}_{\alpha },{\mathcal{H}}_{\alpha }^{\prime }\right);\end{equation} \tag{ 6 }$$

this follows from (3). Furthermore, it can be proven, due to trace preservation of the dynamics, that in every diagonal subspace $\left({\mathcal{H}}_{\alpha },{\mathcal{H}}_{\alpha }\right)$ there is at least 1 stationary state of $\mathcal{L}$ [20]. Hence, the stationary-state degeneracy of a system with a single strong symmetry has a lower bound of n_s.

We now move on to the case where there are multiple, abelian symmetries, i.e. there exists a series of symmetries {S₁, S₂, ..., S_n } whose operator representations R(S_k ) in $\mathcal{B}\left(\mathcal{H}\right)$, which are all linearly independent, satisfy

$$\begin{equation}\left[H,R\left({S}_{k}\right)\right]=0,\quad \left[{L}_{l},R\left({S}_{k}\right)\right]=\left[{L}_{l}^{{\dagger}},R\left({S}_{k}\right)\right]=0\quad \left[R\left({S}_{k}\right),R\left({S}_{{k}^{\prime }}\right)\right]=0,\quad \forall l,k,{k}^{\prime }.\end{equation} \tag{ 7 }$$

Now suppose R(S₁) has n₁ distinct eigenvalues; then we can decompose the Hilbert space into a series of blocks specified by these eigenvalues for this first strong symmetry

$$\begin{equation}\mathcal{H}=\underset{\alpha =1}{\overset{{n}_{1}}{\bigoplus}}{\mathcal{H}}_{\alpha }.\end{equation} \tag{ 8 }$$

We now repeat this procedure, decomposing each H_α into n_α,2 blocks which label the distinct eigenvalues of R(S₂) within the α block. We continue this blocking procedure until we have completely decomposed the space into a series of N spaces. Following our block-decomposition, we then use an identical argument to the previous section and carry this through to $\mathcal{B}\left(\mathcal{H}\right)$, noting the Liouvillian has at least 1 unique steady state in the diagonal spaces.

As a result, the stationary state degeneracy is lower bounded by N, which is the number of permitted distinct combinations of the eigenvalues of the matrices {R(S₁), R(S₂), ..., R(S_n )} (for example there might be eigenvalues of R(S_k ) for which certain eigenvalues of R(S_k') are not possible). We note that the minimum value of N is therefore max({n_k }) and its maximum is ${\prod }_{k=1}^{n}{n}_{k}$ where n_k is the number of distinct eigenvalues of R(S_k ) in the full operator space.

We now derive the main result of this article: the lower bound on the dimension of the nullspace in the case of multiple, non-abelian strong symmetries. We achieve this via the the following theorem and proof.

Theorem 1. Let the dynamics of a system be given by a Lindblad master equation, i.e. equation (1), with the Hamiltonian $H\in \mathcal{B}\left(\mathcal{H}\right)$ and the set of Lindblad jump operators ${L}_{l}\in \mathcal{B}\left(\mathcal{H}\right)$ being bounded linear operators acting over a Hilbert space $\mathcal{H}$. Let there exist a set {S₁, S₂, ..., S_n } of strong symmetries which form a non-abelian group,

$$\begin{equation}\left[H,R\left({S}_{k}\right)\right]=0,\quad \left[{L}_{l},R\left({S}_{k}\right)\right]=\left[{L}_{l}^{{\dagger}},R\left({S}_{k}\right)\right]=0,\quad \forall l,k,\end{equation} \tag{ 9 }$$

where R(S_k ) is the matrix representation of the kth strong symmetry in $\mathcal{B}\left(\mathcal{H}\right)$. In general, due to these symmetries being non-abelian, we have that [R(S_k ), R(S_k')] ≠ 0. We then construct the group representation R via a series of irreducible representations,

$$\begin{equation}R={\bigoplus}_{i=1}^{s}{R}_{i},\end{equation} \tag{ 10 }$$

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 state $\hat{\mathcal{L}}{\rho }_{\infty }=0$ is then bounded from below by ${\sum }_{i=1}^{s}{D}_{i}^{2}$.

This theorem is applicable to both finite and infinite-size systems. However, we note that in an infinite-size system this lower bound is likely to be significantly less than the actual degeneracy due to the unbounded nature of the matrices involved.

Proof. Define ${\tilde {R}}_{i}\left(\cdot \right)$ to be an operator on the full Hilbert space acting as R_i (⋅) in the irreducible representation i and trivially in the the rest of the Hilbert space. Following the decomposition, we have $\left[H,{\tilde {R}}_{i}\left({S}_{k}\right)\right]=\left[{L}_{l},{\tilde {R}}_{i}\left({S}_{k}\right)\right]=0, \forall i,k,l$,. By Schur's representation lemma, this implies that $H,{L}_{l},{L}_{l}^{{\dagger}}$ are multiples of the identity matrix in the R_i block [42]. This, in turn, implies that

$$\begin{equation}\left[H,{\left\{{\tilde {R}}_{i}\left({S}_{k}\right)\right\}}_{nm}\right]=\left[{L}_{l},{\left\{{\tilde {R}}_{i}\left({S}_{k}\right)\right\}}_{nm}\right]=\left[{L}_{l}^{{\dagger}},{\left\{{\tilde {R}}_{i}\left({S}_{k}\right)\right\}}_{nm}\right]=0,\end{equation} \tag{ 11 }$$

where ${\left\{{\tilde {R}}_{i}\left({S}_{k}\right)\right\}}_{nm}$ indicates the matrix at position n and m within the block matrix ${\tilde {R}}_{i}\left({S}_{k}\right)$. This means that there are ${\sum }_{i=1}^{s}{D}_{i}^{2}$ linearly independent operators that commute with $H,{L}_{l},{L}_{l}^{{\dagger}}$. It immediately follows that these operators are in the kernel of ${\hat{\mathcal{L}}}^{{\dagger}}$, which has the same dimension as the kernel of $\hat{\mathcal{L}}$. □

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) can be 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 section 3.2.2. In many situations, however, we find there is no additional structure and that the bound is saturated.

Remark 2. Assuming standard unitary strong symmetries S_k , it is easy to see that [L_l , S_k ] = 0 implies $\left[{L}_{l}^{{\dagger}},{S}_{k}\right]=0$. In this case we can dispense with the requirement $\left[{L}_{l}^{{\dagger}},R\left({S}_{k}\right)\right]=0$ from equation (9) of theorem 1.

Remark 3. Due to the non-abelian nature of the symmetry group there may exist stationary states of $\hat{\mathcal{L}}{\rho }_{\infty ,k}=0$ that are not density matrices, i.e. ${\rho }_{\infty ,k}^{{\dagger}}\ne {\rho }_{\infty ,k}$. However, by nature of the dynamical map in equation (1) every initial density matrix ρ(0) will always remain a valid density matrix: these unphysical states do not corrupt the lower bound we have derived as they are of trace 0 and can always be added to any physical density matrices to make another physical density matrix—provided they preserve its Hermitian nature. The D dimensional kernel always forms a basis for a convex combination of D linearly independent extreme, or fixed, points that are valid density matrices (see [43] for more details). Moreover, these states are interesting from a quantum information and computing perspective because they represent quantum coherences that can be used to store quantum information [44].

In the following, we illustrate these results with a series of examples.

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, 47]. 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 [48]. These types of models have also proven useful in studying excitonic transport in photosynthetic light harvesting systems [49, 50].

Zoom In Zoom Out Reset image size

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

$$\begin{equation}H={\varepsilon }_{\mathrm{g}}\vert 0\rangle \langle 0\vert +\varepsilon \sum _{i=1}^{N}\vert i\rangle \langle i\vert +h\sum _{i\ne j}\vert i\rangle \langle j\vert =\begin{bmatrix}\hfill {\varepsilon }_{\mathrm{g}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \varepsilon \hfill & \hfill h\hfill & \hfill \dots \hfill & \hfill h\hfill \\ \hfill 0\hfill & \hfill h\hfill & \hfill \varepsilon \hfill & \hfill \dots \hfill & \hfill h\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill 0\hfill & \hfill h\hfill & \hfill h\hfill & \hfill \dots \hfill & \hfill \varepsilon \hfill \end{bmatrix},\end{equation} \tag{ 12 }$$

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

$$\begin{equation}\begin{aligned}\hfill {L}_{j}^{+}& ={\mu }_{j}^{+}\vert j\rangle \langle 0\vert ,\hfill \\ \hfill {L}_{j}^{-}& ={\mu }_{j}^{-}\vert 0\rangle \langle j\vert .\hfill \end{aligned}\end{equation} \tag{ 13 }$$

These operators describe the absorption and emission of an exciton on site j where ${\mu }_{j}^{{\pm}}$ 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

$$\begin{equation}{P}_{ij}=I-\vert i\rangle \langle i\vert -\vert j\rangle \langle j\vert +\vert i\rangle \langle j\vert +\vert j\rangle \langle i\vert .\end{equation} \tag{ 14 }$$

The exchange operators that do not act on the bath sites (we label the bath sites as a, b, ....) then form a series of strong symmetries:

$$\begin{equation}\left[H,{P}_{ij}\right]=0,\quad \left[{L}_{a}^{{\pm}},{P}_{ij}\right]=\left[{L}_{b}^{{\pm}},{P}_{ij}\right]=\dots =0,\quad \forall i,j\ne a,b,\dots .\end{equation} \tag{ 15 }$$

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, we can say that every element A formed from different products of exchange elements, and belonging to the permutation group P,

$$\begin{equation}A=\prod _{i,j}{P}_{ij},\end{equation} \tag{ 16 }$$

is a strong symmetry, i.e.,

$$\begin{equation}\left[H,A\right]=0,\quad \left[{L}_{l},A\right]=0,\quad \forall l.\end{equation} \tag{ 17 }$$

It then follows that these group elements are thus left null eigenvectors of the Liouvillian [20, 21]

$$\begin{equation}{\mathcal{L}}^{{\dagger}}\left(A\right)=0.\end{equation} \tag{ 18 }$$

Now suppose we label the sites which are uncoupled form the baths as 1, 2, ..., n. An arbitrary permutation operator can then be block diagonalized as follows

$$\begin{equation}A=\begin{bmatrix}\hfill {R}_{n}\left(A\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {I}_{N-n}\hfill \end{bmatrix},\end{equation} \tag{ 19 }$$

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⟩. This block-diagonal structure should 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 [51]

$$\begin{equation}{R}_{n}\left(A\right)={R}_{\mathrm{t}}\left(A\right)\oplus {R}_{\mathrm{s}}\left(A\right),\end{equation} \tag{ 20 }$$

whose ranks are 1 and n − 1, respectively. Hence, the permutation operator P is further block diagonalized

$$\begin{equation}A=\begin{bmatrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {R}_{\mathrm{s}}\left(P\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {I}_{N-n}\hfill \end{bmatrix}.\end{equation} \tag{ 21 }$$

We have thus found an irreducible representation of the strong symmetries and since the left kernel of $\mathcal{L}$ contains all possible A, by theorem 1, its dimension must be at least (n − 1)² + 1. Furthermore, for every left null vector of $\mathcal{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)² + 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 how the decomposition in equation (20) is reflected in the long-time states of the system. The non-trivial part of A acts on the space spanned by {|1⟩, |2⟩, ..., |n⟩} and to explicitly uncover the decomposition in equation (20) we change from this configuration basis to a new one:

$$\begin{equation}\left\{\vert 1\rangle ,\vert 2\rangle ,\dots ,\;\vert n\rangle \right\}\to \left\{\vert {\phi }_{\mathrm{t}}\rangle ,\vert {\psi }_{1}\rangle ,\dots ,\;\vert {\psi }_{n-1}\rangle \right\},\end{equation} \tag{ 22 }$$

where |ϕ_t⟩ is the trivial representation and is defined as

$$\begin{equation}\vert {\phi }_{\mathrm{t}}\rangle =\frac{1}{\sqrt{n}}\sum _{i=1}^{n}\vert i\rangle ,\end{equation} \tag{ 23 }$$

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 (20) and (21), which provide a natural structure to the steady states of the system.

Zoom In Zoom Out Reset image size

Finally, we note that the Hamiltonian and every Lindblad operators can be likewise block-diagonalized within this decomposition, and take the form:

$$\begin{equation}H=\begin{bmatrix}\hfill \varepsilon +\left(n-1\right)h\hfill & \hfill 0\hfill & \hfill {h}^{{\dagger}}\hfill \\ 0\hfill & \hfill \left(\varepsilon -h\right){I}_{n-1}\hfill & \hfill 0\hfill \\ h\hfill & \hfill 0\hfill & \hfill {H}_{o}\hfill \end{bmatrix},\quad {L}_{i}^{{\pm}}=\begin{bmatrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {L}_{i,o}^{{\pm}}\hfill \end{bmatrix}.\end{equation} \tag{ 24 }$$

Here, we have set, without loss of generality, the ground state energy ɛ_g to zero.

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. Most of the following discussion, however, will be valid for general SU(2) symmetric Hamiltonians. We consider the open isotropic XXX spin chain with the Hamiltonian

$$\begin{equation}H=\sum _{i=1}^{N-1}\left({\sigma }_{i}^{x}{\sigma }_{i+1}^{x}+{\sigma }_{i}^{y}{\sigma }_{i+1}^{y}+{\sigma }_{i}^{z}{\sigma }_{i+1}^{z}\right),\end{equation} \tag{ 25 }$$

where ${\sigma }_{i}^{x,y,z}$ denotes the x, y, z Pauli operator on site i. Importantly, we assume the dissipation in the system commutes with each of the the total spin operators ${S}^{x,y,z}={\sum }_{i=1}^{N}{\sigma }_{i}^{x,y,z}$. If this is the case then the system has an SU(2) symmetry as the total spin operators each form a strong symmetry

$$\begin{equation}\left[H,{S}^{x}\right]=\left[H,{S}^{y}\right]=\left[H,{S}^{z}\right]=0,\end{equation} \tag{ 26 }$$

and are non-abelian as we have [S^α , S^β ] ≠ 0, α ≠ β.

In order to apply our theorem we need to decompose these symmetries over the N-sites into a series of irreducible representations. 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 [52]. We use this notation to perform the required decomposition and as a result 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³)

In order to demonstrate the strength of this bound, we consider coupling explicitly 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 (complex) strength. Such a model should be interesting for understanding effects of dissipative disorder on thermalization and transport properties in a quantum many-body system. The master equatio takes the form

$$\begin{equation}\begin{gathered}{c}\dot {\rho }=-i\left[H,\rho \right]+\sum _{l}\left(2{L}_{l}^{{\dagger}}\rho {L}_{l}-\left\{{L}_{l}^{{\dagger}}{L}_{l}\rho \right\}\right),\\ {L}_{l}=\sum _{i,j}{M}_{l}^{ij}\left({\sigma }_{i}^{x}{\sigma }_{j}^{x}+{\sigma }_{i}^{y}{\sigma }_{j}^{y}+{\sigma }_{i}^{z}{\sigma }_{j}^{z}\right),\end{gathered}\end{equation} \tag{ 27 }$$

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}_{l}^{ij}$ quantifying the strength of this interaction. While 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 the total spin operators form a series of non-abelian strong symmetries.

We test our lower bound by computing the stationary state degeneracy numerically via exact diagonalization. We generate the ${M}_{l}^{ij}$ 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.

Zoom In Zoom Out Reset image size

3.2.2. Total spin dissipation

As a second example of dissipation which preserves the SU(2) symmetry, we consider dissipation induced by the Casimir operator

$$\begin{equation}L={S}^{2}={\left({S}^{x}\right)}^{2}+{\left({S}^{y}\right)}^{2}+{\left({S}^{z}\right)}^{2},\end{equation} \tag{ 28 }$$

such that the master equation reads

$$\begin{equation}\dot {\rho }=-\mathrm{i}\left[H,\rho \right]+\gamma \left(2L\rho {L}^{{\dagger}}-\left\{{L}^{{\dagger}}L,\rho \right\}\right).\end{equation} \tag{ 29 }$$

Since [S², S^x ] = [S², S^y ] = [S², 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.

We compute the actual stationary state degeneracy of this Liouvillian 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. 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² commutes with the Hamiltonian and L is Hermitian there is a basis in which both H and L are mutually diagonal. Clearly, $\hat{\mathcal{L}}$ is also diagonal in this basis. We can write the Casimir operator as,

$$\begin{equation}L=\sum _{i=1}^{s}\sum _{j=1}^{{D}_{i}}\sum _{x=1}^{{M}_{i}}{c}_{i}\vert i,j,x\rangle \langle i,j,x\vert ,\end{equation} \tag{ 30 }$$

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, M₂ = 4 and M₃ = 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,

$$\begin{equation}H=\sum _{i=1}^{s}\sum _{j=1}^{{D}_{i}}\sum _{x=1}^{{M}_{i}}\sum _{y=1}^{{M}_{i}}{h}_{i}\left(x,y\right)\vert i,j,x\rangle \langle i,j,y\vert .\end{equation} \tag{ 31 }$$

We can then diagonalize H in the |i, j⟩ subspace,

$$\begin{equation}H=\sum _{i=1}^{s}\sum _{j=1}^{{D}_{i}}\sum _{x=1}^{{M}_{i}}{E}_{x}^{\left(i\right)}\vert i,j,x\rangle \langle i,j,x\vert ,\end{equation} \tag{ 32 }$$

where ${E}_{x}^{\left(i\right)}$ 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}_{x}^{\left(i\right)}={E}_{{x}^{\prime }}^{\left(i\right)}$.

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}_{x}^{\left(i\right)}={E}_{{x}^{\prime }}^{\left(i\right)}$. Therefore,

$$\begin{equation}\hat{\mathcal{L}}\vert i,j,x\rangle \langle i,{j}^{\prime },{x}^{\prime }\vert =0,\quad \forall x,{x}^{\prime }\vert {E}_{x}^{\left(i\right)}={E}_{{x}^{\prime }}^{\left(i\right)}.\end{equation} \tag{ 33 }$$

Hence, assuming that the energies ${E}_{x}^{\left(i\right)}$ are all distinct the stationary state degeneracy is

$$\begin{equation}{N}_{\mathrm{s}}=\sum _{i=1}^{s}{M}_{i}{\times}{D}_{i}^{2},\end{equation} \tag{ 34 }$$

which is larger than the lower bound of theorem 1. We find that this new prediction for the degeneracy, which accounts for the additional symmetry structure, matches the exact-diagonalization results, see figure 3.

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 chain reads

$$\begin{equation}H=-t\sum _{\langle ij\rangle ,\sigma }\left({c}_{i,\sigma }^{{\dagger}}{c}_{j,\sigma }+H.c.\right)+U\sum _{i=1}^{N}{n}_{i,{\uparrow}}{n}_{i,{\downarrow}},\end{equation} \tag{ 35 }$$

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 [53, 54]. The continuous symmetry of this model is SO(4) and because SO(4) = SU(2) × SU(2)/Z₂ there are two commuting SU(2) symmetries within the model [55].

The first of these is the spin symmetry, generated by the spin operators

$$\begin{equation}{S}^{+}=\sum _{j=1}^{N}{c}_{j,{\uparrow}}^{{\dagger}}{c}_{j,{\downarrow}},\quad {S}^{-}=\sum _{j=1}^{N}{c}_{j,{\downarrow}}^{{\dagger}}{c}_{j,{\uparrow}},\quad {S}^{z}=\frac{1}{2}\sum _{j=1}^{N}\left({n}_{j,{\uparrow}}-{n}_{j,{\downarrow}}\right),\end{equation} \tag{ 36 }$$

which satisfy the SU(2) relations

$$\begin{equation}\left[{S}^{z},{S}^{{\pm}}\right]={\pm}{S}^{{\pm}},\quad \left[{S}^{+},{S}^{-}\right]=2{S}^{z},\quad \left[H,{S}^{z}\right]=\left[H,{S}^{{\pm}}\right]=0.\end{equation} \tag{ 37 }$$

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 [56]. This η-symmetry is generated by the operators

$$\begin{equation}{\eta }^{+}=\sum _{i=1}^{N}{\left(-1\right)}^{i}{c}_{i,{\uparrow}}^{{\dagger}}{c}_{i,{\downarrow}}^{{\dagger}},\quad {\eta }^{-}={\left({\eta }^{+}\right)}^{{\dagger}},\quad {\eta }^{z}=\frac{1}{2}\sum _{i=1}^{N}\left({n}_{i,{\uparrow}}+{n}_{i,{\downarrow}}-1\right),\end{equation} \tag{ 38 }$$

which satisfy the relations

$$\begin{equation}\left[{\eta }^{z},{\eta }^{{\pm}}\right]={\pm}{\eta }^{{\pm}},\quad \left[{\eta }^{+},{\eta }^{-}\right]=2{\eta }^{z},\quad \left[H,{\eta }^{z}\right]=\left[H,{\eta }^{{\pm}}\right]=0.\end{equation} \tag{ 39 }$$

The two SU(2) symmetries are abelian

$$\begin{equation}\left[{S}^{\alpha },{\eta }^{\beta }\right]=0,\quad \alpha ,\beta =+,-,z\end{equation} \tag{ 40 }$$

but clearly the operators which form each SU(2) symmetry are not.

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 describe 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

$$\begin{equation}\left(\mathbf{2},\mathbf{1}\oplus \mathbf{1}\right)\oplus \left(\mathbf{1}\oplus \mathbf{1},\mathbf{2}\right).\end{equation} \tag{ 41 }$$

This means that the four-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

$$\begin{equation}{S}^{x}=\frac{1}{2}\begin{bmatrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{bmatrix},\quad {S}^{z}=\frac{1}{2}\begin{bmatrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{bmatrix},\end{equation} \tag{ 42 }$$

$$\begin{equation}{\eta }^{x}=\frac{1}{2}\begin{bmatrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{bmatrix},\quad {\eta }^{z}=\frac{1}{2}\begin{bmatrix}\hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{bmatrix}.\end{equation} \tag{ 43 }$$

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 (35) to an environment which induces homogeneous, local spin-dephasing. The dynamics is modelled, under the Markov approximation, via the Lindblad master equation

$$\begin{equation*}\dot {\rho }=-\mathrm{i}\left[H,\rho \right]+\gamma \sum _{l=1}^{N}\left(2{L}_{l}\rho {L}_{l}^{{\dagger}}-\left\{{L}_{l}^{{\dagger}}{L}_{l}\rho \right\}\right),\end{equation*}$$

with the Lindblad operators

$$\begin{equation}{L}_{l}={s}_{l}^{z}=\frac{1}{2}\left({n}_{l,{\uparrow}}-{n}_{l,{\downarrow}}\right).\end{equation} \tag{ 44 }$$

A possible experimental implementation for this master equation is discussed in reference [54] 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 while 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

$$\begin{equation}{\mathbf{2}}_{0}\oplus {\mathbf{1}}_{1}\oplus {\mathbf{1}}_{-1},\end{equation} \tag{ 45 }$$

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

$$\begin{align}\hfill & \left({\mathbf{2}}_{0}\oplus {\mathbf{1}}_{1}\oplus {\mathbf{1}}_{-1}\right)\otimes \left({\mathbf{2}}_{0}\oplus {\mathbf{1}}_{1}\oplus {\mathbf{1}}_{-1}\right)\hfill \\ \hfill & \quad \;\;\;\;={\mathbf{3}}_{0}\oplus {\mathbf{1}}_{0}\oplus {\mathbf{2}}_{1}\oplus {\mathbf{2}}_{-1}\oplus {\mathbf{2}}_{1}\oplus {\mathbf{1}}_{2}\oplus {\mathbf{1}}_{0}\oplus {\mathbf{2}}_{-1}\oplus {\mathbf{1}}_{0}\oplus {\mathbf{1}}_{-2},\hfill \end{align} \tag{ 46 }$$

which can immediately be extended to N sites. Using this representation we can apply theorem 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

$$\begin{equation}{N}_{\mathrm{s}}=\sum _{n=1}^{N+1}\left(N+2-n\right){n}^{2}.\end{equation} \tag{ 47 }$$

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⁴), which is faster than the usual O(N³) growth because of the additional U(1) symmetry.