Majorana Scars as Group Singlets

In some quantum many-body systems, the Hilbert space breaks up into a large ergodic sector and a much smaller scar subspace. It has been suggested [1] that the two sectors may be distinguished by their transformation properties under a large group whose rank grows with the system size (it is not a symmetry of the Hamiltonian). The quantum many-body scars are invariant under this group, while all other states are not. Here we apply this idea to lattice systems with N sites that contain M Majorana fermions per site. The Hilbert space may be decomposed under the action of the SO ( N ) × SO ( M ) group, and the scars are the SO ( N ) singlets. For any even M there are two families of scars. One of them, which we call the η states, is symmetric under the group O ( N ) that includes a reﬂection. The other, the ζ states, has the SO ( N ) invariance only. For M = 4 , where our construction reduces to a deformed SU (2) Hubbard chain with local interactions, the former family are the N + 1 η -pairing states, while the latter are the N + 1 states of maximum spin. We generalize this construction to M > 4 . For M = 6 we exhibit explicit formulae for the scar states and use them to calculate the bipartite entanglement entropy analytically. For large N , it grows logarithmically with the region size. In general, the energies of the scars within each family are not equidistant. For M > 6 we also ﬁnd that the scars will local Hamiltonians typically have certain degeneracies.

The past few years have seen growing interest in "quantum many-body scars," the term that was coined in [2]. The essential phenomenon is that there are many-body systems where the Hilbert space breaks up into the bulk of states that satisfy the Eigenstate Thermalization Hypothesis (ETH), and a much smaller scar subspace that does not. Specific constructions of such states have been found in a variety of models . For the recent reviews of the different approaches to scars, see [43][44][45][46].
Remarkably, the quantum many-body scars appear in the commonly used models of condensed matter physics, such as the (deformed) Fermi-Hubbard and t-J-U models on a lattice with N sites. Such models contain two species of complex fermions on each site, c j↑ and c j↓ . In addition to the rotational SU(2) symmetry, they possesses a (broken) pseudospin SU (2) symmetry. The η-pairing states [47,48] form a multiplet of pseudospin N/2, and their role as scars was pointed out and studied in [1,15,21,32]. Another important family are the ζstates that carry the maximum spin N/2; they can be regarded as scars if the SU(2) rotational symmetry is broken [1,32].
There has also been important progress on generalizing the η-pairing states to systems with more than two complex fermions per lattice site [49,50] (for earlier work, see also [51]). In this paper we present a systematic method for constructing multiflavor fermionic systems with weak ergodicity breaking (many-body scars) that leads to generalizing the η and ζ states. It relies on the idea that the scar subspace is invariant under a "large group" whose rank is of the order of the number of lattice sites N [1,32]. In particular, it was noted that the η-pairing states in the spin-1/2 Hubbard model are invariant under SO(N ) which acts on the lattice site index [1], as well as under an even bigger symplectic group Sp(N ) [32]. In this paper we use the idea of group invariance to construct arXiv:2212.11914v1 [cond-mat.str-el] 22 Dec 2022 the scars in lattice systems with an even number M of Majorana fermions per lattice site. This system is equivalent to M/2 complex fermions per site, and for M = 4 we reproduce the previous results singling out N + 1 η-states and N + 1 ζ-states as the scars. For M = 6 we present concise formulae for all the SO(N ) singlets, which come in two families generalizing the η and ζ states. The expressions for the generalized η-states are the same as those found in [50], where it was also understood that they are eigenstates of the SU(3) invariant Hubbard interaction. We use these formulae to calculate the bipartite entanglement entropy analytically. For small N we show numerically that the invariant states we consider possess all the characteristic properties of many-body scars. In particular their entanglement entropy is found to be much smaller than that of the nearby thermal states. For large N we show analytically that the entanglement entropies of scars grows logarithmically with the region size. For M > 4, the scars within the η and ζ families are not in general equidistant in energy even if the Hamiltonian is restricted to local terms only which is a special feature of the present work. We discuss the conditions under which the "revivals" can still be observed.
We also discuss extensions of our approach to M > 6 where the SU(M/2) invariant Hubbard interaction no longer works simply. Therefore, we replace it by another local interaction under which all the SO(N ) singlets are eigenstates. A novel feature we find for the scar states with M > 6 is the presence of degeneracies which typically appear when only local interactions are present, but can be broken by non-local interactions.
We will rely on the method proposed in [1]: the many-body scars span a subspace of the Hilbert space invariant under a large group G that is not a symmetry of the system. The Hamiltonian is chosen to be of the form where H 0 is a simple term governing the dynamics of the scar subspace, T a are generators of the group G, which therefore annihilate the scars, and O a are operators chosen so that the Hamiltonian is Hermitian. In this paper, we apply this construction to the Hilbert space of multi-flavor Majorana fermions on a lattice.

II. DEFORMED HUBBARD MODEL AND SCARS
In this section, we review the SO(4) symmetry of the Hubbard model, using both Dirac and Majorana fermions, and discuss its relation with the scar states in some deformed Hubbard models. For simplicity of the discussion, we consider the model on a 1D lattice of N sites. The standard Hubbard Hamiltonian is the sum of three terms -the hopping, the on-site repulsion, and the chemical potential: Here t is a real hopping parameter, U > 0 is the on-site interaction strength, and c iσ , c † iσ are the fermionic ladder operators satisfying the anticommutation relations The magnetic field is (µ ↑ −µ ↓ )/2 while the standard chemical potential is (µ ↑ + µ ↓ )/2. We find it convenient to perform a site dependent phase rotation, c jσ → e −ij π 2 c jσ , upon which the hopping term acquires an imaginary coefficient: In these variables, the Hubbard Hamiltonian is H Hub = T + V + µ.
The time reversal symmetry in this case is defined as follows For the vanishing magnetic field, µ ↑ = µ ↓ , the Hamiltonian H Hub has both time reversal symmetry and the spin SU(2) symmetry, which acts on the spin index σ. The generators of SU(2) are where ψ A j , A = 1, 2, 3, 4 are four Majorana fermions on site j. Then H Hub at µ ↑ = µ ↓ = U 2 admits a manifestly SO(4) invariant form where the hopping terms A ψ A j ψ A j+1 can be identified as special cases of the (antihermitian) SO(N ) generators T jk ≡ where the η-vacuum |0 η is the same as the empty vacuum |0 , and the ζ states The η states span an (N + 1) dimensional representation of the pseudospin group SU(2) and the ζ states furnish a spin N 2 representation of the spin group SU (2). The ζ states have fixed fermion number N and have eigenenergies −(µ ↑ − µ ↓ )m − µ ↓ N with respect to the Hubbard Hamiltonian H Hub . The η states are also energy eigenstates of the Hubbard Hamiltonian, i.e. H Hub |m η = m(U − µ ↓ − µ ↑ )|m η . In [32], the authors constructed the Hubbard model deformed by quartic OT terms that break both the spin and pseudospin symmetries. Then the SO(N ) invariant η and ζ states remain eigenstates, and they have all the typical properties of scar states. In the next section we will consider a different set of symmetry breaking deformations which also keep the η and ζ states as scars. We will also extend the construction from 4 Majorana fermions per site to a higher even number M .

III. MULTI-FLAVOR MAJORANA FERMIONS ON A LATTICE
Generalizing the Majorana description of the Hubbard model, we consider a lattice of N sites, hosting an even number M flavors of Majorana fermions ψ A j , A = 1, 2, · · · , M on each site. Their anti-commutation relations are invariant under the action of SO(N ) × SO(M ) group. We can build (antihermitian) generators of SO(N ) and SO(M ) out of these fermions Their commutation relations are given by Grouping the Majorana fermions into complex combinations provides a convenient way to construct states in the Hilbert space. On each site j, define α = 1, . . . , M/2 flavors of complex fermions which satisfy the standard anticommutation relations {c iα , c † jβ } = δ αβ δ ij . Following the general recipe described in appendix A, we construct a basis of so(M ) in terms of these complex fermions that makes its root decomposition structure manifest:  or their direct sums. When N is even, it was found in [54] that the SO(N ) singlets furnish the reducible representation λ + N/2 ⊕ λ − N/2 by using a character method. We will present an elementary way to show that this structure holds for both odd and even N in section IV, where we study more detailed structures of these singlets.
A. Ergodic Hamiltonians that support singlet states as many-body scars Following the recipe proposed in [1], we construct (local) Hamiltonians of the form H = H 0 + O ij T ij . The first term H 0 is designed to leave the space of SO(N ) singlets (denoted by S) invariant. The second term, which will be referred to as the OT term henceforth, should be hermitian and break some symmetries of H 0 , by choosing O ij properly. The Hamiltonians we discuss below should have as few symmetries as possible such that they are guaranteed to produce ergodic bulk spectra and such that the many-body scars are not fully occupying an isolated sector of a symmetry. Nevertheless since we consider bosonic Hamiltonians the fermion number parity symmetry is always preserved The hopping term T M is summed over nearest neighbors. It can be thought as a special OT term with O ij ≡ it when In complex fermion basis, H M becomes where n jα = c † jα c jα is the fermion number operator of flavor α at site j. Since n jα takes value in {0, 1}, we have 1 − 2n jα = (−1) njα and hence V M can be simplified as where n j is the total fermion number at site j. Turning on chemical potentials for each flavor, e.g.
breaks the SO(M ) symmetry of H M . On a general ground, one can consider the Hamiltonian H µ = AB iµ AB J AB , where µ AB is a real anti-symmetric matrix. But one can always diagonalize the matrix µ AB and redefine the Majorana fermions to get (III.12). However, the SO(M ) Casimir C

SO(M ) 2
is still a conserved charge of the combined Hamilto- In the Hubbard model, generic chemical potentials µ ↑ and µ ↓ break the SO(4) symmetry. In this case, the Casimir operators of SU(2) and SU(2) are conserved charges because the symmetry breaking term is a linear combination of η 3 and ζ 3 . The OT term can be used to break the conservation of C . For instance, the following sextic interacting term does this job where r AB are randomly chosen real numbers. This form of local interacting Hamiltonian is analogous to that adopted by Shiraishi and Mori [4], j P j h j P j , where P j is a set of local projection operators that satisfy P 2 j = P j . While our SO(N ) generators T jk are not projectors and do not commute, our singlet conditions on the scar states, T jk |φ = 0 are analogous to the conditions P j |φ = 0 imposed in [4]. In the recent papers [39,40], some parallels were drawn between the approach identifying scars as the sector invariant under a large group [1,32] and the Shiraishi-Mori construction.
Adding up H M , H µ and H int , we get the following Hamiltonian and it produces a good separation between scars and other states in the entanglement entropy plots. Alternatively, instead of the sextic interaction term (III.13), we may consider the quartic term with O jk being a hermitian and satisfying [O jk , T jk ] = 0. It is easy to check such H OT is hermitian and annihilates the singlets. There is a simple way of construction of such operators O jk . For that purpose let us note that T jk could be regarded as generator of a U (1) Q T symmetry. This gives the following charges An operator O jk that has zero charge Q T then commutes with the hopping operator T jk . For instance, we can consider the following term, that is biliniear in d jk,A Another natural generalization of the Hubbard model would be on-site density-density interaction between different flavors, which can be described by the following potential [55] This term does not have SO(M ) symmetry for M ≥ 6. It has been noticed in [50] that although V M still keeps S invariant when M = 6, it can typically map SO(N ) singlets to nonsinglets when M ≥ 8.

B. Controlling the position of scars in the spectrum
One of the possible strategies [1] to control the position of scars in the spectrum relies on the addition of a term H C that annihilates all the scars but acts positive-definitely on all other states. Because scars are SO(N )-invariant in our case the most obvious choice for such a term is the quadratic Casimir operator of the SO(N ) group where T ij are SO(N ) generators given in eq. (III.2). The interaction in eq. (III.19) is however highly non-local in real space. We find numerically that for the accessible system sizes also the following much more local interaction can be used where the summation is only over the nearest neighbours It is non-negative definite, and because of the presence of SO(N ) generators, we have where |s is an SO(N ) singlet state and |ns is a non-singlet state. Therefore, it can be used to move the scars in energy with respect to all other states without changing the relative position (and the revivals period) of scars themselves. In particular, using eq. (III.20) one can achieve that the low-energy part of the spectrum is comprised of many-body scars only.

IV. SO(N ) SINGLETS AS SCARS
In this section we discuss various properties of the SO(N ) singlets. In particular, for M = 6, i.e. 6 Majorana fermions per site, we can write down the wavefunction of every SO(N ) singlet. Similar explicit formulae for the η states have been obtained in the SU(3) Hubbard model [50]. Furthermore, we show analytically that both the η states and the ζ states have a sub-volume law for entanglement entropy in the thermodynamic limit.
Denote the irreducible representation containing |0 by H η , which corresponds to the highest-weight vector λ η = N 2 , N 2 , · · · , (−1) M/2 N 2 , and denote the irreducible representation containing |0 ζ by H ζ , which corresponds to the highest-weight vector repeatedly on the empty vacuum |0 , and hence are generalizations of η states in the Hubbard model. Similarly, H ζ is built with The two representations H η and H ζ are always distinguished by a reflection operator in O(N ), which can be realized by e iπnj . H η has parity +1 under e iπnj and hence is O(N ) invariant. H ζ , on the other hand, has parity −1 under e iπnj and hence is only SO(N ) invariant. The fully occupied state, which is apparently a highest-weight state, has parity (−1) M/2 under the reflection e iπnj . So it belongs to H η when M/2 is even and H ζ when M/2 is odd. This also explains why the highest weight vectors of H η and H ζ depend on the parity of M/2. When N is odd, H η and H ζ can also be distinguished by the fermionic parity (III.8). In this case, H η has fermionic parity +1 and H ζ has fermionic parity −1. Altogether The spectrum of the chemical potential term H µ restricted to S is encoded in the partition function Z(β) ≡ tr S e −βHµ , which group theoretically is equal to the sum of two SO(M ) characters where x α ≡ e −βµα . Each character can be computed using Weyl character formula (A.8). These Weyl characters admit the following expansions where {p α } take values in all integers when N is even and all half integers when N is odd. The coefficients D η p (D ζ p ) are nonnegative integers, and they vanish for all but a finite number of vectors p. Because of the identification x α = e −βµα , the characters implies that there are D η p (D ζ p ) linearly independent states in H η (H ζ ) that diagonalize all {h α } simultaneously with eigenvalues {p α }. Altogether, the eigenenergies of the full Hamiltonian restricted to the scar subspace can be summarized as For generic chemical potentials µ α and interaction strength U , the energy spacing does not have a common divisor. So we will not observe revivals starting from a generic state in S. However, because {p α } are integers or half-integers, revivals are possible with special choices of µ α . For example, the revival can appear when We conjecture that the remaining degeneracies we observe for M > 6 are "unbreakable," i.e. they cannot be removed by any local perturbations that preserve the decoupling of the scars. On the other hand, it is not hard to see that they can be broken by the non-local interactions, such as where r AB are a set of real random numbers.

C. Product scar states
Product states are very special because their entanglement entropy vanishes. In S, product states are either empty or fully filled for each of the M/2 flavors. Hence there are 2 M/2 such states. To describe their wavefunctions, we define the following N -fermion operators Then the 2 M/2 product states can be expressed as where 1 ≤ α 1 < α 2 < · · · < α κ ≤ M/2. In particular, κ = 0 corresponds to the empty vacuum |0 , and κ = M/2 corresponds to the fully occupied state. Since A † α has odd parity under the reflection e iπnj , the product state |α 1 , · · · , α κ belongs to H η when κ is even and H ζ when κ is odd. For the full Hamiltonian H, c.f. eq. (III.14), |α 1 , · · · , α κ is an eigenstate with energy In the case of M = 4, the product states in S are where the κ = 0 and κ = 2 states are η states |0 η and |N η (c.f. eq. (II.9)), and the two κ = 1 states are ζ states |0 ζ and |N ζ (c.f. eq. (II.10)), which have total spin ± N 2 .

D. The M = 6 case
In the Majorana model of six flavors, we have a very good analytical control of the singlet states, from eigenstate wavefunctions to entanglement entropy.

Wavefunctions and eigenenergies
When M = 6, the O(N ) invariant subspace H η has dimension N +3 3 . Consider states in H η of the form (N −k12−k13−k23)! and hence is nonvanishing when k T ≡ k 12 + k 13 + k 23 ≤ N . Next, these states are linearly independent, because {k 12 , k 13 , k 23 } uniquely fixes the fermion numbers of the three flavors, namely k 12 + k 13 particles of flavor 1, k 12 + k 23 particles of flavor 2 and k 13 + k 23 particles of flavor 3. Finally, counting nonnegative integer solutions of the inequality k 12 + k 13 + k 23 ≤ N precisely reproduces N +3 3 . Therefore, a normalized orthonormal basis of H η that diagonalizes the three Cartan generators {h 1 , h 2 , h 3 } simultaneously is where the normalization factor C k (N ) is given by These states (IV.14) are also constructed in [50] as energy eigenvectors of the 3-component Hubbard model. The eigenenergy of |k 12 , k 13 , k 23 with respect to the Hamiltonian H is 16) For the SO(N ) invariant subspace H ζ , we can construct an orthonormal basis similarly They are eigenstates of the Hamiltonian H with energy 18) where µ 1,2 = µ 1,2 , µ 3 = −µ 3 .

Analytic computation of the entropy of the singlet states
With explicit wavefunctions known for the SO(N ) singlets in the M = 6 case, we proceed to compute their entanglement entropy analytically. We will focus on H η but the same recipe of computation also works for H ζ . We want to mention that a similar calculation is also done for η states in usual Hubbard model in [57], and for a special class of η states in multicomponent Hubbard model in [50].
We divide the underlying lattice into two disjoint subsets Σ 1 and Σ 2 . For example, Σ 1 consists of the first N 1 sites, i.e. i = 1, 2, · · · , N 1 , and Σ 2 consists of the rest N 2 = N − N 1 sites. In each sublattice Σ a , there is an empty vacuum |0 a satisfying |0 = |0 1 ⊗ |0 2 . Then we are allowed to split the η-operators in the following way such that η a † αβ can excite η states of the subsystem Σ a on |0 a Because of η † αβ = η 1 † αβ + η 2 † αβ , we have a tensor product decomposition of any η state |k ≡ |k 12 , k 13 , k 23 defined in eq. (IV.14) let us note that |m 1 and |k − m 2 are also singlet states in each subsector. Taking the partial trace over the Hilbert space of Σ 2 yields the reduced density matrix ρ Σ1 (k) of |k The density matrix ρ Σ1 corresponds to a pure state if (i) all k αβ are vanishing, or (ii) one k αβ is equal to N and the rest are vanishing. The former is trivial since it implies m 12 = m 13 = m 23 = 0. For the latter, say k 12 = N , we have first m 13 = m 23 = 0. Then nonvanishing λ k (m) requires m T ≥ k T − N 2 = N 1 and m T ≤ N 1 , which completely fix m 12 = N 1 . Indeed, (i) corresponds to |0 , and (ii) corresponds to A † α A † β |0 , 1 ≤ α < β ≤ 3, which are the only product states in H η .
Next we proceed to compute the entanglement entropy of ρ Σ1 (k) in the thermodynamic limit, defined as the limit of N → ∞, k αβ → ∞ such that ν αβ ≡ k αβ N are finite. Since k T < N , ν αβ should satisfy ν T ≡ ν 12 + ν 13 + ν 23 < 1. We further choose N 1 N so that Σ 2 can be treated as a heat bath and meanwhile keep N 1 1 to allow scaling of entanglement entropy ρ Σ1 . In this limit, λ k (m) is sharply peaked around m αβ = m * αβ ≡ ν αβ N 1 . At this extremal , and for general m, where M is a 3 × 3 symmetric matrix, given by with determinant equal to ν 12 ν 13 ν 23 (1 − ν T ). The scaling property of the entanglement entropy of ρ Σ1 (k) in thermodynamic limit can then be computed by replacing the sum over m with a triple integral d 3 m: This integral shows that S Σ1 (k) scales as the logarithm of the size of Σ 1 . The coefficient 3 2 for the logarithmic behavior agrees with [50].
The whole derivation of entanglement entropy above also works for the other SO(N ) invariant subspace H ζ . In particular, starting with a ζ state |k ζ , we end up with same density matrix (IV.20), with |m 1 being replaced by the corresponding ζ states on Σ 1 .

E. The structure of singlets for M ≥ 8
When M = 8, there are six different η † αβ . In order to construct explicit wavefunctions of H η , we consider the following set of states, which generalizes the M = 6 case, , which is smaller than N +3 Noticing that |0 is actually the highest-weight state of H η , we build upon it another set of states which has dimension N +5 6 . Since states in V I η have total fermion number Q ≤ 2N and states in V II η have total fermion number Q > 2N , the two sets V I η and V II η have no intersection. Adding up their dimensions gives exactly the dimension of H η . Therefore an orthonormal basis of H η is where k αβ and αβ satisfy α<β k αβ ≤ N, α<β αβ < N .
Similarly, one can show that H ζ is an orthonormal direct sum of V I ζ which is spanned by We will compute their entanglement entropy numerically in the next section. For M ≥ 10, it becomes very hard to write down explicit wavefunctions of all SO(N ) singlets. In [50], some special states belonging to H η are considered, namely (η † 12 ) k2 · · · (η † 1,M/2 ) k M/2 |0 with k 2 + · · · + k M/2 ≤ N . They are eigenstates of the generalized Hubbard potential (III.18) and their entanglement entropy can be evaluated analytically. In the thermodynamic limit, the entanglement entropy is found to scale as S Σ1 ∼ M −2 4 log N 1 [50].

F. Off-Diagonal Long-Range Order
By construction the singlet states are quite non-local and we can check that by studying the spatial dependence of the following operators It is easy to check that if i = j this operator is real and does not depend on the indices i and j. Indeed, let us consider a rotation Q ik that acts in the (i, k) plane and leaves the rest vectors untouched. Then it is easy to check that the value of this operator depends on the choice of AB and the singlet |s ∈ S. But if we sum over A and B we would get a simpler operator where we have used anticommutation relations. After that we can get where we have used the fact that |s is annihilated by the action of hopping T ij = A ψ A i ψ A j . Another "sum rule" arises if we average O AB ij over all the singlet states which amounts to computing the trace of ψ A i ψ B i ψ B j ψ A j over the scar subspace. When i = j, its value is independent of the choice of i, j. For A = B, we have which implies that O AB (s) is nonvanishing for at least one |s .

V. NUMERICAL RESULTS
To test the above predictions numerically we implement the Hamiltonian where the auxiliary symmetry-breaking term H SB is given by either  As the first test we examine the dimension of the SO(N )singlet subspace. To do this we scale the part (T M + H SB ) of the Hamiltonian (V.1) that is proportional to SO(N ) generators. Then we count the number of energy levels that remain unchanged up to numerical precision. The number of these states is given in the Table V.1 for all accessible system sizes and agrees with the number of SO(N ) singlets given in eq. (IV.3).
Many-body scars are expected to stand out by violating the eigenstate thermalisation hypothesis in that an observable measured in these states largely deviates from the thermal average at the same energy or temperature. We use entanglement entropy as the observable of interest. To calculate it we imagine cutting the system in the middle such that all the states left from the cut can be associated with system Σ 1 and right from the cut with system Σ 2 . We perform the Schmidt decomposition (numerically performed as a singular value decomposition) of an eigenvector with respect to this separation: where the states |x j Σ1 and |x j Σ2 describe only the left/right sub-systems respectively and Λ j are singular values such that Λ 2 j are the eigenvalues of the reduced density matrix of either subsystem (their analytical expressions are known for M = 6 eq. (IV.21)).
The entanglement entropy is then calculated as and quantifies the amount of entanglement between the subsystems. The minimum value of 0 is achieved when the full state |x is a product state.
In addition we will also study the statistics of the level spacings in the spectrum. Going through the full spectrum we determine the level spacing and level spacing ratio for every pair of energy levels.
Mean values for this level spacing ratio is known analytically [58] for certain types of random matrices: r ≈ 0.5359 for the generalized orthogonal ensemble (GOE, real matrices) and r ≈ 0.6027 for the generalized unitary ensemble (GUE, complex matrices).

A. M = 4
The M = 4 case discussed in Sec. II has the same Hilbert space as the spin-1/2 electron models studied in Ref. [1] and features two families of scars. The entanglement entropy in the full Hilbert space and the level spacing histogram in the even sector of the fermion number parity (III.8) are shown in Fig. V.1 for a system with N = 8, M = 4 and open boundary conditions. Eq. (III.13) was used as a symmetry-breaking term.
Both families of scars have the entanglement entropy significantly lower than the generic states at the same temperature. Two states (n = 0 and n = N ) in each scars family are product states and have exactly zero entanglement entropy.
Level spacings are analysed in the even sector of the fermionic number parity, they are qualitatively the same in the odd sector. The average ratio of r ≈ 0.60038 we obtain is very close to the GUE value. The histogram of level spacings (excluding the singlet energy levels) shown in the bottom panel of Fig. V.1 shows that the near-zero gaps are almost absent in the spectrum as is expected due to level repulsion and absence of residual symmetries. Altogether, we observe chaotic bulk spectrum and the singlet states clearly violating ETH which therefore attests that the singlet states are manybody scars in this system. The average level spacing ratio in the even sector of the fermion number parity is 0.5998 and is approximately the same as for hermitian random matrices (GUE). Vanishing probability of near-zero gaps that can be seen in the bottom panel of Fig. V.2 indicates level repulsion characteristic of ergodic quantum systems without symmetries.
Combining these observations we can conclude that also in the M > 4 case the SO(N ) singlets have all the properties of the many-body scars.
We note that with the choice of "random" µ α the system in question has scars (their energy given in eq. (IV.7)) that are not equally spaced in energy. Such situations have not been described in literature to our knowledge.
The off-diagonal long-range order (eq. (IV.31)) measured between the most distant sites (i = 1, j = 4) for A = 1 and B = 2 is shown in Fig. V.3. As predicted by the eq. (IV.35) the sum normalized by the dimension of the scar subspace equals 1 28 and is independent of the choice of A and B. We confirm numerically that the full minimal set of ODLRO measurements contains M − 1 operators (IV.31) where A and B are chosen as (k, k + 1) for k between 1 and M − 1. This means that at least one of the M − 1 correlators is non-zero in every scar state.

VI. DISCUSSION
We have presented the structure of the quantum manybody scars in 1D lattice systems of N sites with M Majorana fermions per site. Following the idea of group-invariant scars [1,32] we identified the classes of Hamiltonians where the SO(N ) invariant states are exact eigenstates and showed that they have the properties expected of the scars. We obtained a variety of analytical and numerical results, including the explicit wave functions for M = 6, which for the generalized η-pairing states were found also in [50]. We used the wave functions to calculate the entanglement entropies analytically; as in the earlier work [50,57], they exhibit logarithmic growth with the system size.
The Hilbert space of multi-flavor Majorana fermions used in this work allowed us to uncover several possibilities for the behavior of scars that were not discussed in the earlier  literature. The number of states breaking ergodicity in a system with N sites and M flavors grows as N M (M −2)/8 . For M > 4 this is faster than the linear growth that was a feature of the scars discussed previously. Furthermore, for M > 6 we find degeneracies in the scar subspace that cannot be lifted by the local interactions which preserve the decoupling of scars. These degeneracies present a new promising resource that could potentially be used for robust quantum computing, similarly to how the topological degeneracies are used in topological quantum computing schemes [59].
Although the many-body scars we study can exhibit revivals under some conditions, they are not in general equallyspaced in energy. This is the first such example which departs from the concept of a large spin precessing in a magnetic field that accounted, to our knowledge, for all the examples of the scars to date except for the ones constructed using the Shiraishi-Mori method [4].
Finally, the multi-Majorana systems possess a rich variety of the off-diagonal long-range order correlators. The result of measuring such correlators in a scar state does not depend on the distance between the points. By analogy with the superconducting correlations found in the η states for M = 4, the full set of the corresponding ODLRO operators we identified for M > 4, may find multiple interpretations and applications in further studies. where the overall normalization factor 1 2 is inserted such that [ζ † αβ , η † βγ ] = η † αγ . The generators corresponding to negative roots −(e α ± e β ) are hermitian conjugate of ζ † αβ and η † αβ . Altogether, the root decomposition of so(2n, C) is so(2n, C) = h 1≤α<β≤n Cζ † αβ ⊕ Cη † αβ ⊕ Cζ αβ ⊕ Cη αβ (A. 5) An integral highest-weight vector λ can be parameterized by λ = n α=1 λ α e α , where λ α are either integers or halfintegers, satisfying λ 1 ≥ λ 2 ≥ · · · ≥ λ n−1 ≥ |λ n |. In terms of Young diagram, λ α is the number of boxes in the α-th row. Given a highest-weight vector λ = (λ 1 , · · · , λ n ), the Casimir C