Group theoretic approach to many-body scar states in fermionic lattice models

It has been shown [arXiv:2007.00845] that three families of highly symmetric states are many-body scars for any spin-1/2 fermionic Hamiltonian of the form $H_0+OT$, where $T$ is a generator of an appropriate Lie group. One of these families consists of the well-known $\eta$-pairing states. In addition to having the usual properties of scars, these families of states are insensitive to electromagnetic noise and have advantages for storing and processing quantum information. In this paper we show that a number of well-known coupling terms, such as the Hubbard and the Heisenberg interactions, and the Hamiltonians containing them, are of the required form and support these states as scars without fine-tuning. The explicit $H_0+OT$ decomposition for a number of most commonly used models, including topological ones, is provided. To facilitate possible experimental implementations, we discuss the conditions for the low-energy subspace of these models to be comprised solely of scars. Further, we write down all the generators $T$ that can be used as building blocks for designing new models with scars, most interestingly including the spin-orbit coupled hopping and superconducting pairing terms. We expand this framework to the non-Hermitian open systems and demonstrate that for them the scar subspace continues to undergo coherent time evolution and exhibit the"revivals". A full numerical study of an extended 2D $tJU$ model explicitly illustrates the novel properties of the invariant scars and supports our findings.


I. INTRODUCTION
In recent literature there has been considerable interest in the many-body scar states  (for their pedagogical overviews, see [38,39]). The scar states are typically spread throughout the energy range of the spectrum and are, therefore, relevant when the excitation energy is high. They do not obey the eigenstate thermalization hypothesis (ETH) [40][41][42], one of the most fundamental conjectures that allows us to bridge quantum mechanics with statistical physics. On the practical side, an initial state made of many-body scar states repeatedly returns to itself in time evolution preventing the loss of quantum information through thermalization. This offers intriguing prospects for quantum computing.
In our previous paper [1] (see also [23,24]), we presented a general strategy for systematically designing the Hamiltonians with a many-body scar subspace S invariant under the action of a continuous group G, which is bigger than the symmetry group of the Hamiltonian. The general form of such Hamiltonians is H = H 0 + j O j T G j , where T G j are generators of the symmetry group G and O j are a set of operators such that the product O j T G j is Hermitian. H 0 must admit the states in S as eigenstates and the revivals are observed when the gaps between the corresponding eigenvalues have a common divisor.
In this work, we demonstrate that many of the commonly used condensed matter models of interacting electrons actually happen to be of this form and specify their group-invariant scar subspaces S. This applies to the Hubbard, Heisenberg, and some other interactions, and any models constructed out of them on various lattices and in arbitrary dimension such as the extended 2D tJU model that we consider in detail as a prototype. Most of these models can be readily implemented in experiment paving the way for their systematic experimental studies. Our results show that some of the known scar states [15,21] can indeed be explained by the mechanism of Ref. [1]. An important feature of our construction is that G is a large Lie group, such as O(N ) or SU(N ), whose rank is of order N , the number of lattice sites. The generators of such groups include the nearest neighbor hopping terms, and using them we can construct a variety of (approximately) local lattice models.
We also extend the framework of [1] to non-Hermitian Hamiltonians. In particular, we study a model where the product OT is not Hermitian. This opens the way to studying the coherent time evolution of a group-invariant scar subspace in open systems. To our knowledge, many-body scars at finite energy density in non-Hermitian systems have not been discussed previously [49], while the non-stationary, periodic phenomena in dissipative strongly-interacting systems is an emergent hot topic [50][51][52][53][54][55].
In Sec. II we specify the three families of scar states that are relevant to our discussion of spin-1 2 fermionic models. We describe their structures, and the generators T annihilating them. These generators may be viewed as the free Hamiltonians. In Sec. III we show how the well-known interaction terms, such as the Hubbard and Heisenberg interactions, decompose in terms of the same generators. As a consequence, the models including any linear combinations of such interaction terms also have H 0 + OT decomposition. In Table III we give these decompositions for several models, such as the J 1 − J 2 and Haldane-Hubbard models. We conclude with an example (Sec. V ) where the full construction is detailed for  the tJU model and the numerical evidence of many-body scar states is provided.

II. GROUP GENERATORS AND INVARIANT STATES
Consider a lattice of N sites with each one hosting a complex spin-1/2 fermion (electron). Their two spin states are created on site j by operators c † jσ , where σ =↑, ↓. We do not assume any spatial structure in general but will specify it where necessary. It is useful to think of the Hilbert space of 2N fermionic degrees of freedom as factorized according to the representations of one of the following groups:  [56,57]. Due to the peculiar nature of the Hilbert space, the representations of a group acting on the site index are locked to particular representations of the dual O(4) = SU(2) spin × SU(2) η [56] (the SU(2) η group is often called the pseudospin). The symmetry properties of the invariant states we are going to consider in this work are summarised in Table I (the states are given explicitly in (22), (24), (25)). Note that each of these subspaces is invariant under two groups, a large group acting on the site index and a small group acting on the spin or pseudospin. For example, the states |n ζ are invariant under both U(N ) and SU(2) η . Therefore, the generators of each of these groups T U(N ) and T SU(2)η and the commutators [T U(N ) , T SU(2)η ] all annihilate [58] the |n ζ states and can be used as T in the H 0 + OT Hamiltonian. Further insights into the interrelations between the various groups and the invariant states are provided in Fig. 1 and Appendix C.
In order to construct the Hamiltonians such that the invariant states are many-body scars, we now discuss the generators of all the relevant groups. We will use the standard notation for the fermion numbers at site i: The total fermion number is and we will call the corresponding symmetry U(1) Q (the actual generator is (Q − N )/2). The generators of the rotation group SU(2) spin are given by The spin operator at site i is where σ A are the Pauli matrices, and the Greek indices take two values, ↑ and ↓. In particular, is the magnetization at site i. The symmetry corresponding to total magnetization is U(1) M ⊂ SU(2) spin . Another important group is the pseudospin [59][60][61], which is denoted by SU(2) η . Its generators are The generators of SU(2) spin and SU(2) η involve summation over all lattice sites. Therefore, using such generators in the interaction terms j O j T G j will produce a very non-local Hamiltonian. Instead, we will follow the suggestion in our previous paper [1] to construct these interaction terms using the generators of large rank groups, such as SU(N ), which include the nearest-neighbor spin-independent hopping terms.
The full set of generators of SU(N ) are the hopping terms with generally complex amplitudes: There is an SO(N ) subgroup of SU(N ) whose generators are the spin-preserving hopping terms with imaginary amplitudes They are invariant under SU(2) spin × SU(2) η .
An alternate basis for the generators of SU(N ) is where t a are the traceless Hermitian N × N matrices. The simple root generators of this algebra could be chosen as the nearest-neighbor hoppings Interrelations between the three invariant subspaces and the corresponding groups. The three groups that act naturally on the Hilbert space are shown inside the ellipses, the corresponding generators are shown inside the rectangular boxes, and the invariant states are shown inside the circles.
where n(i) is a nearest neighbour of the site i on a given lattice. Commuting simple roots we can restore the whole algebra (8). Any Hermitian spin-independent hopping terms on any lattice are linear combinations of T ij (8). They belong to the SU(N ) algebra and annihilate the SU(N ) singlets (and consequently also the singlets of the larger group U(N )).
There is another, less obvious, group SU(N ), whose generators are (9) and the spin-orbit coupled hopping The spin-preserving hopping terms (9) are generators of the SO(N ) which is a subgroup of both SU(N ) and SU(N ) (see Fig. 1). By combining operatorsT kl (12) and Q A (4) we obtain more general spin-orbit coupled hopping terms that annihilate the singlets of U (N ) and SU(2) spin (|n η states). A special case of these generators is the local spin operator S A i . By commuting η + (7) and T ij (8) we can get the operators that annihilate the singlets of SU(2) η and of U(N ) (|n ζ states) There is a particular set of generators of U(N ) that we will find useful: 1  2  3  4  5  6  7  8   9  10  11  12  13  14  15  16   17  18  19  20  21  22  23  24   25  26  27  28  29  30  31  32   1  2  3  4  5  6  7  8   9  10  11  12  13  14  15  16   17  18  19  20  21  22  23  24   25  26  27  28  29  30  31  32 FIG. 2. A depiction of the bipartite square 2D lattice. The red lines indicate the "even" hoppings that connect sites on different sublattices. The green lines indicate the "odd" hoppings that connect sites on the same sublattice. The numbers labelling the sites show the lattice site index, which is acted on by one of the large groups such as SU(N ).

For example,
where T N 2 −1 is the generator of SU(N ) corresponding to the matrix We have N i=1 K i = 2η 3 , and the square of K i can be expressed in terms of the local magnetization: where we used n 2 i↑ = n i↑ and n 2 i↓ = n i↓ . Now let us consider a special case of a bipartite lattice where the vertices are divided into two non-intersecting sets, which we can call red and blue (see Fig. 2). Here, in addition to the above operators, we can define the transformed genera-torsT where and the sum runs over the blue sites. The generators (19) form a group that we denote as SU(N ) . A particular subset of the generatorsT is the nearest-neighbor hopping terms with a real coefficient Other ways of constructing a SU(N ) algebra of longer-range hopping terms are discussed in the next sub-section. We will denote the spin preserving subgroup of this group as SO(N ) ⊂ SU(N ) . We now provide the basis for each invariant subspace listed in the Table I. Each subspace should be invariant under the action of H 0 for the invariant states to be scars [1]. A particular H 0 may have eigenstates that are different from the basis discussed below but may be obtained from it by a rotation.
The subspace invariant under U(N ) and SU(2) η is spanned by |n ζ which have the highest possible physical spin; namely, they form the spin-N 2 representation of SU(2) spin [1]: where n = 0, . . . , N , and is the spin raising operator. The states |n η are invariant under U(N ) and SU(2) spin ; they form the N +1 dimensional representation of pseudospin SU(2) η : where η = η + (7). On a bipartite lattice we can further define the states |n η that are invariant under U(N ) and form the N + 1 dimensional representation of SU(2) η : These states are known as the η-pairing states [59,62] [63]. It is not hard to check that Combined with the fact that both families of states include the vacuum, that is annihilated by S A i , this means that all the |n η and |n η states are annihilated by the local spin operators S A i . The importance of the spin and pseudopspin SU(2) groups has been previously understood in the context of the Hubbard model [59][60][61]. These two SU(2) groups are related by the Shiba transformation, that acts on the fermion operators in the following way: We note that this transformation interchanges the |n η and |n ζ states and leaves the SO(N ) generators (9) invariant. The Hubbard Hamiltonian is also known [64] to be invariant under the Shiba transformation.  All three families of scars have a very high degree of symmetry with respect to the spatial index j. As a consequence, when the sites are arranged into a lattice, all three families are translation and permutation-invariant. In Table II we summarize the action of all the generators we mentioned on the three families of the invariant states. This table provides a simple way of constructing Hamiltonians with one of the families as scars -any linear combination of the operators with 0 in a corresponding row can be used as an operator T in the general H 0 + OT form. In Sec. III we show that many of the commonly known interaction terms, such as the Hubbard and Heisenberg interactions, can be decomposed in terms of the generators listed above and therefore have invariant manybody scar states.
A. Longer-range hopping terms Some classic models, e.g. the Hubbard model, include only the nearest-neighbor hopping term with a real coefficient t: However, the models that aim to describe realistic materials, such as the high-T c superconductors, may include the nextto-nearest-neighbor or even longer-range hopping terms with real coefficients. It is, therefore, interesting to inquire if the deformation by such additional hopping terms preserves the families of scar states. In this section we discuss when such terms are generators of the relevant symmetry groups and annihilate the paired states |n η . Any Hermitian spin-independent hopping terms on any lattice (including T , T O ij ) are special cases of the complexamplitude hopping T ij (8). Therefore, they belong to an SU(N ) algebra and annihilate the |n ζ subspace. The nearestneighbor hopping with real amplitude (21) on a bipartite lattice is a generator of SU(N ) and annihilates |n η .
A complication with the longer-range real-amplitude hopping comes from the fact that a commutator of such terms may not always be expressed as T r kl alone. Yet, in some special cases, a constraint on the lattice type, hopping range or boundary conditions allows this algebra to be closed. For example T ij is restricted to the nearest-neighbor hopping on a bipartite lattice where we can perform a transformation (20) on c k which maps t → it. Then the hopping terms become generators (9) of SO(N ), which is important for the argument that the states |n η are scars.
However, if the next-to-nearest-neighbor (nnn) hopping terms are added, the transformation keeps those hopping coefficient real. As a result, the O(N ) invariant states |n η are not scars in presence of the nnn hopping terms. Indeed, combining the nearest-and next-nearest-neighbor hopping does not lead to a closed algebra. For example, let us consider the following set of hopping terms connecting sites 1, 2, 3 in a linear chain.
The hoppings T 1 and T 2 are nearest-neighbor, while T 3 is next-nearest-neighbor. Then It is not of the form of T 3 ; therefore, the subalgebra of hopping terms with real coefficients does not close.
To avoid this issue, we may use a bipartite lattice and restrict the allowed hopping terms. A closed SO(N ) algebra is formed by the "even" long-range real-amplitude hopping that contains exclusively the terms connecting the sites belonging to different groups of the bipartite lattice (red-to-blue) in Fig.  2, meaning that it is a hopping over 2k neighbours, where k is a non-negative integer. The number of the nearest neighbors hopped over by the term t ij σ c † iσ c jσ is determined by the lattice as the smallest number of connected sites that need to be visited when travelling from site i to site j. The evennearest-neighbour hopping is shown in red lines in Fig. 2 and includes the nearest-neighbour hopping as a subset. Such hopping annihilates |n η as a generator of SO(N ) and |n ζ as a generator of SU(N ).
The algebra SO(N ) can also be constructed on a bipartite lattice when both even-and odd-nearest-neighbour hopping are present in the system. In this case, we must require that all the even-nearest-neighbor hopping terms have real amplitudes while all the odd-nearest-neighbor terms have imaginary amplitudes (for example, in Fig. 2 we may include the nearestneighbor hoppings with a real amplitude and the shortest diagonal hoppings with an imaginary amplitude). Such a hopping will again annihilate both |n ζ and |n η subspaces. An example of this scenario is the Haldane model [65] discussed in Sec. A 6.

B. Electromagnetic field
In the absence of an electromagnetic field, the hopping terms are typically taken to have real coefficients. The presence of a (possibly time-dependent) electromagnetic field in a lattice model introduces the phase factors e iαij c † iσ c jσ . As discussed above, such spin-independent hopping terms with a complex amplitude T ij are generators of SU(N ) and therefore annihilate the |n ζ states.
The only effect of an electromagnetic field on the |n ζ states comes from the coupling between the magnetic field and the spin of the electrons: δH = −µ Q · B, which results in a linear in B and equidistant splitting of the |n ζ states according to the projection of their total spin and changes the "revivals" period. Other families of the invariant states are mixed by the electromagnetic field with the rest of the spectrum. Because |n ζ forms the maximum-spin representation of SU(2) spin these states get split the most by the magnetic field w.r.t. to all other states in the Hilbert space. If all other terms in the Hamiltonian are bounded this means that in a strong enough magnetic field an |n ζ state can always be made the ground state.
For some particular configurations (for example a closed 1D chain with a magnetic field of π through it), the modified hopping amplitude may be purely imaginary. In such cases, the hopping T O ij (9) is a generator of SO(N ), which makes the states |n η insensitive to the magnetic field (they also have zero spin).

III. INTERACTION TERMS
Let us consider an arbitrary lattice with N sites on which spin-1 2 electrons are placed, and let i and j refer to two sites on this lattice. We would like to rewrite the commonly used electron-electron interaction terms in the form H 0 + OT and analyze their action on the three invariant subspaces.

A. Hubbard interaction
In the Hubbard model, two electrons interact only when they are located on the same site i. The Hubbard interaction at site i may be written in terms of the generators K i of group U(N ): This annihilates the |n ζ states. While neither |n η nor |n η are eigenstates of (32), they are eigenstates of the Hubbard interaction summed over the full lattice that we write using (32), (18) and (15) as Since the local magnetization M i annihilates the states |n η and |n η , we have which implies that the energies of these scar states have equal spacing. This is related to the fact that the pseudospin generators η ± do not commute with the Hamiltonian of the Hubbard model [59][60][61].
A slightly modified Hubbard interactioñ can be cast in the form and has the advantage that all three invariant families are eigenstates on every site:H Hub

B. Density-density interaction
The density-density interaction at different sites, H dd ij = n i n j , can also be written in terms of the generators (15): For the SU(N )-invariant states |n ζ , we have H dd ij |n ζ = |n ζ . Neither |n η nor |n η states are eigenstates of H dd ij since K i does not act on them with a definite value.
The generalized density-density interaction reads and M i |n η = M i |n η = 0 (see Tab II). M i does not annihilate |n ζ and mixes them with the non-singlet states. However |n ζ are eigenstates of the total magnetization: Similarly, |n η and |n η are not eigenstates of n i , but for the total fermion number Q we have Q |n η = 2n |n η and Q |n η = 2n |n η . For |n ζ we have n i |n ζ = |n ζ . The above observations mean that H gdd ij annihilates |n η and |n η states when α ij = 0 and leaves |n ζ unchanged when β ij = γ ij = δ ij = 0; for example α ij n i n j |n ζ = α ij |n ζ .

C. Heisenberg interaction
The Heisenberg interaction H Heis ij = S i · S j couples spins on two different sites, i and j. Using the expression (5) for the spin operator in the Hilbert space of spin-1/2 complex fermions, we find This may be written as where and we introduced the off-diagonal generators of SU(N ): Eq. (40) represents the H 0 +OT decomposition of the Heisenberg interaction w.r.t. to the group U(N ). It follows that the U(N )-invariant states |n ζ are degenerate eigenstates of S i · S j on any lattice, in any dimension and for any i = j with energy 1 4 . We also note that the Heisenberg interaction annihilates the states |n η and |n η because they are annihilated by the local spin operators S A i (see Eq. (26) and the discussion following it).

D. Symmetry-breaking perturbation
To highlight the ergodicity-breaking properties of the invariant states and to be able to tune to the fully chaotic regime we will consider a simple symmetry-breaking term of the OT form that we write down for a rectangular lattice in two dimensions, where i labels the horizontal and j the vertical direction: where r ij , q ij ∈ [0, 1] are real random numbers and where r M , q M are also real random numbers. Both S hor ij and S vert ij are special cases of the hopping T ij (21) which is simultaneously a generator for SU(N ) and SU(N ) groups (see Sec II and Tab. II). Therefore the full perturbation term H p H is of the pure OT form and annihilates two invariant families: |n ζ and |n η . It can be added to any model supporting these states as scars without changing their energy but breaking all symmetries except the particle number conservation.
We note that the symmetry-breaking terms such as (43) also lead to the absence of any dynamical symmetries that could otherwise cause persistent oscillations of local observables [66].

E. Non-Hermitian perturbation
A non-Hermitian Hamiltonian may result under certain approximations from reducing a closed full system to an effective description of an open sub-system. The non-Hermitian Hamiltonian may not conserve the norm of the state which corresponds to the probability leaking out or into the open system. While the derivation of such an effective description is beyond the scope of this work we show here that the invariant states remain stable and decoupled also in non-Hermitian systems as long as the non-Hermitian Hamiltonian has the form H 0 +OT . We consider the following non-Hermitian OT term This only differs by a minus sign from Eq. (43), and thus also annihilates the same two invariant families |n ζ and |n η . This operator is invariant under complex conjugation, and its eigenvalues form complex conjugate pairs. According to the classification of Ref [67,68], the model falls into the Ginibre symmetry class AI, and therefore will have the same level statistics as the Ginibre GUE (GinUE) distribution [69]. Another possibility to obtain a non-Hermitian term of the OT form annihilating invariant states is by taking only a "half" of the hoppings (8), (9), (12), (13), (21), without the hermitian conjugate. A 2D modified Hubbard model of this form where t 1 = t 2 are real parameters, was recently shown [70] to be amenable to quantum Monte-Carlo simulations. Since such non-Hermitian hopping terms are certain generators of the SU(N ) group, they annihilate the |n ζ states.

IV. MODELS WITH INVARIANT SCARS
The interaction terms considered above, as well as their linear combinations, are of the form H 0 +OT . Together with the generators from Sec II they can be used as building blocks for designing Hamiltonians in which some of the invariant states |n ζ , |n η , and |n η are many-body scars.
Many commonly used models, such as the Hubbard, Heisenberg, and tJU models fall into this class. In Table III we explicitly re-write them as H 0 + OT and indicate the invariant states comprising the scar subspace, the derivations and some more details are given in Appendix A. The presence of |n η scars in the extended Hubbard and Hirsch models, that has been demonstrated in the literature [15,21], is a direct consequence of the group theoretic structure presented in this work.
Because several families of scars with different symmetries may be present in the same model simultaneously the H 0 + OT decomposition is different between |n ζ and the paired |n η / |n η states. However, the difference is only in that a certain constant term (shown in third column of Table  III) is either a part of H 0 or OT depending on which group we consider (full separate expressions with respect to each group are given in Appendix A).
To facilitate the design of custom models not listed here we also added Table IV to the Appendix which gives the action of other common Hamiltonian terms on the invariant states.
Note that in models that conserve spin, the |n ζ states that comprise a maximum spin-N/2 representation of SU(2) spin will not be true scars. To make them such, adding a perturbation breaking spin-conservation (such as (43) or (45)) is required. Analogously, the states |n η only become true scars upon addition of an η-pairing symmetry breaking perturbation of the OT form.
A. Engineering the location of scars in the spectrum

Making a scar the ground state
The energies of the scars and the basis in the invariant subspace are determined by the H 0 part of the Hamiltonian, and we can change their position in the spectrum by adding a term to H 0 that commutes with it (is diagonal in the basis selected by H 0 ).
The |n ζ states form the maximum spin representation of SU(2) spin and have a definite total spin and its axis projection quantum numbers. This means that a sufficiently strong magnetic field will make the |n ζ state with the largest axis projection the ground state. It also controls the splitting between the |n ζ states and the revivals period. For the basis in Eq. (22) one would use the B z = Q 3 (4) magnetic field while Haldane-Shastry Well-known models written as H0 + OT . The decomposition is different for paired |n η /|n η and |n ζ states. The term in the third (H ζ 0 = OT η ) column is part of H0 for the decomposition with respect to the U(N ) group (|n ζ scars) and a part of OT with respect to U(N ) or U (N ) groups (|n η or |n η scars). Eqs. (18) and (40) allow to easily switch between the expressions for |n ζ and paired states. All models are considered in the unrestricted Hilbert space that allows any on-site occupations incl. double-occupations. All decompositions are valid in any dimension and on any lattice except for the models that involve the |n η states -they only become scars when the lattice is bipartite. N nn 1 and N nn 2 are the numbers of nearest-and next-nearest-neighbour pairs in a particular lattice. Cij (41) are generators of SU(N ). The last three columns indicate which families of the invariant states are many-body states in the corresponding model. *: only when Φ = π/2, t1 = 0. **: only when Φ = π/2 for the basis (54) used later in our numerical example we will use the B y = Q 2 field (4).
The states |n η and |n η have a definite particle number and include the states with maximal and minimal possible particle number (all filled |N η or |N η and all empty |0 η or |0 η ). This guarantees that a sufficiently large chemical potential term can be used to make one of these two states the ground state.

Coupling to the rest of the scar subspace
Suppose we are at zero temperature and are in the ground state which per the above results is a scar state. Now we turn on for a limited time the "raising operator" of the corresponding scar family.
In case the ground state is an all-filled |N η or |N η state we add a term η − + η + to the Hamiltonian (see (25) and (7)). This may potentially be made by placing the system into direct contact with a superconductor. As a result, the system will be initialized to a state that is a linear combination of |n η or |n η . Moreover, since |n η form an irreducible representation of SU(2) η we can by the action of a group e it(αη − +α * η + ) send any initial state, say |0 η or |N η , to any desired linear combination of the paired states.
In case the ground state is a state with maximal spin |N ζ we add a term ζ + ζ † (23) to the Hamiltonian. This may be implemented by introducing an external magnetic field (B x = Q 1 for basis (22) and B z = Q 3 for basis (54)). As a result the system will be initialized to a state within the |n ζ subspace. The same algebraic argument as for |n η above shows that with an appropriate choice of magnetic field and interaction time one can turn any initial state in the |n ζ subspace into an arbitrary desired linear combination of |n ζ states.
Finally, we can mix the states from the families |n η and |n ζ by the simultaneous action of magnetic field and η term.
After the above steps the system is initialized to a state in the singlet scar subspace. At this point one may choose to lower or turn off the magnetic field or the chemical potential that was used in the first step IV A 1 to make one of the scar states the ground state.
It would be interesting to analyze the relation of this scheme to a number of protocols for preparing a paired |n η state that rely on dissipation or periodical driving [73][74][75][76].

Low-energy subspace composed of scars only
We can also arrange for all the low-energy subspace be composed solely of many-body scars. To do that we need to add a non-negative definite operator P to the Hamiltonian that annihilates the scars (and the scars only). Increasing the magnitude of such a term will leave the scars untouched but will push the rest of the spectrum up in energy.
For the scar states invariant under a particular group we have various options for designing the desired operator P . For example, P can be chosen to be the quadratic Casimir operator of the corresponding group, as was done in [77]. However, from the point of view of lattice models this is a complicated and non-local operator. A simple choice of a local operator is P = k (T k T † k ) l , where T k are the simple root generators of the group (11) and l > 0 is an integer.
The Hubbard interaction already includes (33) a similar term P = i K 2 i with K i the generator of U(N ) (15). It does annihilate the |n ζ scars and pushes non-scar states high in energy for large U at half-filling (where the H 0 energy contribution of the Hubbard interaction is zero). However it also "accidentally" annihilates some of the non-scar states and thus doesn't alone allow to create a scar-only low-energy subspace.

V. TWO-DIMENSIONAL tJU MODEL
We illustrate the concepts discussed above using the example of a perturbed tJU model on a 2D rectangular bipartite lattice shown in Fig. 2. The Hamiltonian of the standard tJU model [80][81][82]

combines the Hubbard and Heisenberg interactions
and can be viewed as a generalization of Hubbard or t − J models relevant for high-T c superconductivity [83]. It is typically assumed that t is real and negative, so that the kinetic energy is minimized at zero momentum. Certain types of longer-range hopping could be considered in addition without changing the structure of the scar subspace as discussed in Sec. II A.
The H 0 + OT decomposition of the model (47) with real t can be performed w.r.t. to two groups U(N ) and U(N ) which leads to two families of group-invariant scars |n ζ and |n η . In case of U(N ) we have while for the group U(N ) it is  (48) and (49) are only different by assigning a certain constant term to H 0 or OT as indicated in the decomposition given in Table III. We recall that the states invariant under group G become scars in a model that can be written as H 0 + k O k T k , where T k are generators of G. Note that a generator of any subgroup of G is also a generator of G and can also appear as T k in the decomposition. This is actually the case in (49) (and all other decompositions w.r.t. paired states) where the hopping terms are generator of SU(N ) , while the Heisenberg and magnetization terms involve the local spin S i that is also a part of the full symmetry group of |n η (see Sec. C).
The standard tJU model (47) conserves the total physical spin. Therefore the states |n ζ form a separate symmetry sector of this model. To make them true scars we break the total spin conservation by adding a perturbation (43). The full Hermitian Hamiltonian we study numerically reads where we added a term proportional to the SU(2) spin generator Q 2 (4). It acts as H 0 on |n ζ and splits them according to the index n: Q 2 |n ζ = (2n − N ) |n ζ ). For the SU(2) spininvariant states |n η it is of the pure OT form and annihilates them exactly. Note that by increasing γ we can make the scar state with maximum Q 2 , |S 1 (54), the ground state. This may be used to initialize the system to a state from many-body scar subspace as described in Sec. IV A. The full non-Hermitian Hamiltonian we consider is where for numerical investigations we set β 1 = 0.4β. In both the Hermitian and the non-Hermitian cases, the part of the full Hamiltonian that acts with a constant on the invariant states is The states |n η are already eigenstates of (53) as written in Eq. (25) while the states in Eq. (22) are not the eigenstates of (52). Instead, the basis in the SU(N )-invariant subspace determined by the Hamiltonian (52) reads The energies of the invariant states are given by  where n -is the index of a state in its respective family (25) or (54).
The energy of the product state |0ζ tJU = |S 1 in large systems is proportional to N and therefore has good chances to be the ground state at half-filling for (J/2−µ−γ) < 0. An |n η state is more likely to be the ground state when (U − 2µ) < 0.
Because both Hermitian and non-Hermitian models we consider have an exact H 0 + OT decomposition, they have the two families |n ζ and |n η as scars for any choices of the coupling constants. The two scar families form two equidistant towers of states with the energies given in (55). Revivals within each individual subspace can be observed for any values of the couplings. However, to see the revivals of an initial state that is a mix of |n ζ and |n η subspaces all the gaps between them must have a common divisor which represents a constraint on the choice of the constants µ, U , and γ.

A. Numerical results
For the numerical experiment we use the rectangular 3×3 lattice and choose open boundary conditions to make sure the real-amplitude nearest-neighbor hopping is a generator of U(N ) and annihilates the |n η states (boundary condition wouldn't matter in a system with even number of sites). We set t = −0.4, J = 1, U = 8, µ = 1.3 and r M = 1.426974 and q M = 2.890703 (see Eq. (44)). For γ = 1 this corresponds to the g.s. filling of ν = Q 2N = 0.44(4), 11% below the half-filling which corresponds to the potentially high-T c -relevant regime [84,85]. For our simulation we instead choose γ = 3.6. At this value the half-filled |S 1 state (54) becomes the ground state. This simplifies the initialization of the system to the scar subspace in experiment, which is discussed in detail in Sec. IV A. Because |S 1 is a product state, it can alternatively be created by application of a simple gate circuit on each site (see SM of Ref [1]).

Quantum chaos
The level statistics parameters of the hermitian model r = 0.5306 (r GOE = 0.5359) and the non-hermitian model r nh = 0.7378 (r Ginibre ≈ 0.74) are close to the values of the corresponding random ensembles (defined in Appendix B) and thus indicate that the bulk spectra of both systems are fully chaotic. This is further elaborated by the gap distribution P (s) and by the presence of the "dip-ramp-plateau" structure in the spectral form factor (Fig. 3 and 4 ) typical for chaotic systems.
The level spacings distribution, P (s), is the probability density function of the spacings between consecutive eigenvalues. A key feature of the P (s) of random Gaussian matrices is that it decays to zero as s → 0. This phenomenon is called level repulsion and implies that it is very unlikely for eigenvalues to be nearly identical. We observe level repulsion in both Hermitian (Fig. 3) and non-Hermitian (Fig. 4) systems studied.
The spectral form factor SFF is usually defined as and gives information about the longer-range correlations of eigenvalues. The main elements of the SFF for a random matrix is a dip ramp plateau structure. The ramp is caused by the repulsion of eigenvalues that are far apart; these eigenvalues are anti-correlated, which is why the ramp is below the plateau. The plateau is a result of generic level repulsion, as degeneracies are unlikely. The ramp and plateau occur at later times and thus probe shorter distances, and these elements are a result of a phenomena known as spectral rigidity. The dip occurs at early times and so it probes larger distances; it is the Fourier transform of the entire spectrum. The probability density function (pdf) of the eigenvalues for a random Hermitian matrix follows a semicircular distribution. This shared property of the pdf leads to a similar slope of the dip in the SFF of the Wigner random ensembles. In the Hermitian case the dip ramp plateau structure is seen at high enough temperature (Fig. 3). In the non-Hermitian system the dip ramp plateau structure is best seen if the SFF is calculated for the real parts of the eigenvalues only. In both systems the large magnetic field γ = 3.6 (used to make one of the scar states the ground state) causes the correlations at one corresponding frequency that results in the peak seen in the SFF plot soon after the dip. This peak is absent or much less pronounced for moderate magnetic field of γ = 1 (see Appendix B).

Spectrum
In the spectrum of the non-Hermitian Hamiltonian (51) shown in Fig. 5 we observe that all the scar states remain at the real axis and are not effected by the non-Hermitian terms while the majority of the eigenvalues of the non-Hermitian Hamiltonian become complex. This demonstrates the stability of the many-body scar states in suitably designed open systems. The non-Hermitian spectrum has a "conjugation symmetry": for every state with energy a + ib there is another state with energy a − ib. All the observables measured in any two such states (such as entanglement entropy) are also equal. For this reason we choose to plot such observables as a function of the real part of the energy eigenvalue: Re E.

Entanglement entropy
The entanglement entropy at half-filling is shown in Fig.  6. We observe that in both Hermitian and non-Hermitian systems the |n ζ scar states expectedly have entanglement entropy much lower than the rest of the spectrum at corresponding energy (temperature). This is also true for the |n η scars that are located in the respective particle number sectors (not shown). High magnetic field used to make one of the |n ζ states the ground state leads to the emergent structure most clearly seen in Fig. 6 for the non-Hermitian case: each total spin sector acquires different energy and starts to form a separate arc as it would be if the total spin was exactly conserved by the Hamiltonian. This additional structure disappears for smaller magnetic field of γ = 1 as seen in the Appendix Fig.  13.

ETH violation
To further demonstrate the violation of strong ETH we evaluate the "superconducting" and "magnetic" off-diagonal long-range order (ODLRO) correlators characteristic of the |n η [59] and |n ζ [1] states respectively. We observe (see Fig. 7 and 8) that the corresponding expectation values are significantly different in the invariant states relative to the micro-canonical average in both systems which allows us to conclude the invariant states are scars in this system.
Because of the high symmetry of the invariant states, the above correlators, when evaluated for scars, do not depend on the coordinates x 1 , y 1 , x 2 , y 2 [1]. For the numerical evaluation we set the points 1 and 2 to be the most distant points in our system with open boundaries: (x 1 , y 1 ) = (1, 1) and (x 2 , y 2 ) = (3,3).
Note that while all the |n ζ states are at half-filling (Q = 9) the particle number of the paired states |n η is 2n and therefore only one such state is visible for any fixed filling (the data shown is for Q = 8).
Very strong magnetic field is present in both systems and couples to the states with non-zero magnetization. This results in the spikes seen in the data for the non-Hermitian system which is apparently more susceptible to the magnetic field.

Time evolution and revivals
One of the most striking and counter-intuitive features of the scar states in the non-Hermitian Hamiltonian (51) is the stable and coherent time evolution of the scar subspace shown in Fig. 9. The system is initialised to a state that is a uniform mix of all the 2N + 2 = 20 scar states present in the system. As shown in Fig. 10 this state is coming back to itself exactly after the time intervals 2π/ω ≈ 20.94, where ω = 0.3 is the greatest common divisor of all the gaps between the scar states according to (55). The norm of the state is preserved throughout, although the system is open (Hamiltonian is non-Hermitian).
For the initial state composed of generic eigenstates the imaginary components of the eigenvalues lead to the probability density quickly flowing into the effectively open sys-tem. As can be seen in the left panel of Fig. 10, the norm of the time-evolved state and correspondingly the overlap explode already over one revivals period.
A similar phenomenon is observed also if only 1 percent (by weight) of generic, complex eigenstates is admixed to the initial state dominated (99 percent by weight) by scars. The initially small imaginary components obtain exponential amplification with time, however a few periods of the revivals with the same period as in other systems can still be observed as shown in the right panel of Fig. 10.
As expected, starting from a mix of scar states, stable revivals are also observed for the Hermitian Hamiltonian (see Fig. 11). In contrast, the information initially stored in the generic states quickly dissipates through thermalisation (right panel of Fig. 11).

VI. DISCUSSION
We have shown explicitly that the Hamiltonians of some well-known models (Table III) are of the H 0 + OT form; therefore, they support the group-invariant scars for any coupling constants and without a need for fine-tuning. A large number of other models can be built by combining the terms we list as the group generators in Sec. II. Of particular interest are the "superconducting" terms in (14) and the spin-orbit coupled hopping terms in (12) and (13). It will be very interesting to explore the interplay of the weak ergodicity breaking and the superconductivity or topology by studying the models built of such terms.
Another possible generalization of our approach is to consider discrete groups G instead of Lie algebras. In this case we should use T = g − 1 with H = H O + OT , where g ∈ G and g −1 H 0 g = H 0 .
From the quantum computation perspective, the groupinvariant scar states we consider are analogous to a "decoherence-free" subspace that can be used to reduce noise in the universal quantum computation [86,87]. Thus, another interesting direction for future work is the development of specific protocols that would enable more robust quantum computation.
Many of the models we are considering are among the most commonly and widely used in condensed matter physics. Therefore there is a great potential for studying these models and the weak ergodicity breaking effects therein in cold atoms quantum simulations, and also for identifying materials that could host this phenomenon. We do describe conditions under which the ground state is a many-body scar or the full low-energy subspace consists solely of scars, which should further facilitate the experimental explorations.
While in this work we considered the case of spin-1/2 fermionic systems, the general approach of Ref. Hamiltonian. This greatly expands the realm of weak ergodicity breaking phenomena and will hopefully inspire new theoretical and experimental studies. In particular, it would be interesting to study the relation of the degenerate groupinvariant scar subspaces in the non-Hermitian systems (e.g. in (55), B y = 0 for |n ζ or µ = U/2 for |n η ) to the "exceptional points," the degeneracies in the non-Hermitian operators that lead to numerous exotic phenomena and are being intensively studied currently [54,[88][89][90]. are the numbers of nearest-and next-to-nearest-neighbor pairs in a particular lattice.

Hubbard model
The Hubbard model Hamiltonian reads where i and j label sites of an arbitrary lattice in arbitrary dimension, i, j stands for the nearest neighbors, and t is a negative real number. Substituting our result (32) we obtain the H 0 + OT form in terms of the generators of U(N ), where the hopping term is a generator of SU(N ) as a special case of T ij (8). On a bipartite lattice this hopping term coincides with (21) is therefore a generator of SU(N ) . Using (18) we obtain the H 0 + OT decomposition for the group SU(N ) : Comparing the two expressions (A1) and (A2), we notice that they are similar except that the constant U ·N/2 belongs to H 0 in the first and to OT in the second equation (as a consequence of (18)).
Only the |n ζ states are scars for the Hubbard model on an arbitrary lattice. On a bipartite lattice, because the H 0 + OT decompositions are possible w.r.t. two different groups SU(N ) and SU(N ) , the singlets of these groups |n ζ and |n η together form the 2N + 2-dimensional scar subspace.
The states |n η are not eigenstates of the Hubbard model; they are mixed by the hopping terms.

Heisenberg model
The spin-1/2 Heisenberg Hamiltonian is The constraint S i · S i = 3/4 translates into the requirement of half-filling on each site: n i = 1. This sector of Hilbert space contains the |n ζ states and a single, half-filled |n η / |n η state. When we restrict to the locally half-filled subspace K i = 0, we get from (40) and the H 0 + OT form (w.r.t. SU(N )) is where N nn 1 is the number of nearest-neighbour links on a given lattice and E ij are the SU(N ) generators.
Therefore the U(N )-invariant states |n ζ are degenerate eigenstates of the Heisenberg model with the energy E Heis
We have seen that each of such terms acts on |n ζ states with a constant (Sec. III C and Eq. (40)): S i · S j |n ζ = 1 4 |n ζ . Therefore all the |n ζ states are degenerate in the isotropic Haldane-Shastry model, are scars and have the energy which is always integer in units of π 2 4N 2 . We note that Ref. [91] does mention the existence of such integer-energy states. [93] Both families of paired states are exactly annihilated by the Haldane-Shastry Hamiltonian as a special case of the Heisenberg interaction. 4. J1 − J2 model J 1 -J 2 model is the Heisenberg interaction between nearestneighbours (J 1 ) and next-nearest-neighbours (J 2 ): Using Eq. (40) we write where the SU(N )-generators C ij are defined in Eq. (41). It follows that H J1J2 |n ζ = 1 4 (J 1 N nn 1 + J 2 N nn 2 ) |n ζ , while H J1J2 annihilates the paired states |n η and |n η just as each of its terms does individually. Therefore all three families decouple in the J 1 -J 2 model and are scars.
We also note that the above results directly apply to the Majumdar-Ghosh model which is a special case with J 1 = 2J 2 .

Hirsch model
a. The "reduced" version [15,71] The "reduced" Hirsch model for which the |n η states were shown to be scars in Ref. [15] reads where X is a real number and −σ is the opposite direction of spin σ or equivalently [15] H Hir Defining noticing that the terms T ij are generators of both SU(N ) and SU(N ) (see (21)) and using Eq. (32) for n j,↑ n j,↓ we obtain the H 0 + OT form for the model (A14) with respect to the group SU(N ) (|n ζ scars) and with respect to the group SU(N ) (|n η scars) On a generic lattice T ij is a generator of SU(N ) and the half-filled SU(N )-invariant |n ζ states become scars with the energies, H Hir r |n ζ = −µN |n ζ . (A18) On a bipartite lattice the nearest-neighbour hopping T ij is a generator of SU(N ) . Therefore in addition to the |n ζ the |n η states become scars with the energies H Hir The full version from Ref. [72] The Hirsch model in its original formulation (Eq. 6 in Ref. [72], see also [94]) contains an additional density-density term but is missing the chemical potential Defining and using (37) for density-density we bring the full Hirsch model to the H 0 + OT form w.r.t. the group SU(N ) We have seen in Sec. III B that the generic density-density interaction can be written in terms of the SU(N ) generators but doesn't have the paired states as eigenstates.
Therefore only |n ζ states with energies H Hir |n ζ = V N nn 1 |n ζ (A23) are scars in the full Hirsch model [72].

Haldane-Hubbard model
The spin-full Haldane-Hubbard model is [65,95,96] H HH = t 1 i,j ,σ c † jσ c iσ + c † iσ c jσ − µQ+ and can be defined on any bipartite lattice and in particular on the 2D honeycomb lattice of graphene. The model is obtained by adding the complex-amplitude next-nearest-neigbour hopping to the Hubbard model which means the H 0 +OT decomposition only differs from Hubbard model by the addition of the t 2 hopping to the OT (see Tab. III).
In the most general case the t 2 hopping terms are the generators of SU(N ) and the model possesses only one scar family |n ζ with the corresponding H 0 + OT decomposition For Φ kl = π/2 the amplitude of the t 2 hopping becomes purely imaginary and because the lattice in the Haldane-Hubbard model is always bipartite the hopping t 1 + t 2 is a generator of SU(N ) (see Sec II A) and the |n η family of scars is added while the corresponding H 0 + OT decomposition reads If in addition t 1 = 0 then the remaining t 2 hopping alone is a generator (see Eq. (9)) of SO(N ) which is a subgroup of both SU(N ) and SU(N ) and the third scar family of SU(N ) × SU(2) spin -invariant states |n η is added. The corresponding H 0 + OT decomposition is

Appendix B: Quantum chaos
We can better understand several properties of our system using knowledge from random matrix theory. There are three main random matrix ensembles corresponding to the Hermitian models that we study in this paper; these are commonly known as the Wigner ensembles consisting of the Gaussian Orthogonal Ensemble (GOE), the Gaussian Unitary Ensemble (GUE), and the Gaussian Symplectic Ensemble (GSE). The GOE is time reversal invariant and is a random real symmetric matrix (H = H T ) where the entries are drawn from a normal Gaussian distribution. The GUE is not time reversal invariant and is a random Hermitian matrix (H = H † ) where the entries are drawn from a complex Gaussian distribution. Finally, the GSE is time reversal invariant (but breaks rotational symmetry) and is comprised of real quaternion matrices. Here, we quantify quantum chaos in our models using three distinct measures, and we can compare some of these measures from our model to those of the corresponding random matrix ensemble. We examine the spectral form factor (SFF) [97], the level spacings distribution (P (s)), and the mean spacings ratio, r .

a. Hermitian Hamiltonian
The mean level spacings ratio, r , is often used to quantify chaos as well as spectral transitions between Wigner-Dyson ensembles. The spacing ratio, r, is defined as, where λ i is the i th eigenvalue. See equation (B2) for the definition for a non-Hermitian system. We also see level repulsion in the probability density function of the consecutive level spacings ratio, P (r), as we send r → 0. The analytic mean level spacings ratios are calculated in [78], and are r ≈ 0.5359 , r ≈ 0.6027, and r ≈ 0.6762 for the GOE, GUE, and GSE respectively. Based on the Figs. 3 and 12, we can conclude that our Hermitian models have a chaotic bulk; the SFF has a dip ramp plateau structure, and the level spacings plots and r values closely match those of the GOE. We note that the peak is absent or much reduced for the moderate magnetic field γ = 1 (see Fig. 12).

b. Non-Hermitian Hamiltonian
The Ginibre symmetry classes are the non-Hermitian analogs to the Dyson symmetry classes. We can compute the level spacings of our non-Hermitian models and compare to those of the Ginibre random matrix analog. It is also possible to compute complex spacing ratios r . For example, the Ginibre GUE (GinUE) r is numerically determined to be ≈ 0.74 [79].
The r-value is defined as [79].
where E n,1 and E n,2 are the nearest and the next nearest energy levels in the spectrum to the given energy level E n . We can see from Fig. 4, that our non-Hermitian model shows evidence of a chaotic bulk. The level spacings plot fits closely to that of the Ginibre distribution, and the r value of our model is also close to the Ginibre value. The interpretation of the SFF for the non-Hermitian model is less straightforward, though we do see a dip ramp plateau structure when considering only the real part of the eigenvalues. The dip ramp plateau structure is less clear when considering only the magnitude or imaginary part of the eigenvalues.