A monomial matrix formalism to describe quantum many-body states

We propose a framework to describe and simulate a class of many-body quantum states. We do so by considering joint eigenspaces of sets of monomial unitary matrices, called here ‘M-spaces’; a unitary matrix is monomial if precisely one entry per row and column is nonzero. We show that M-spaces encompass various important state families, such as all Pauli stabilizer states and codes, the Affleck–Kennedy–Lieb–Tasaki (AKLT) model, Kitaev's (Abelian and non-Abelian) anyon models, group coset states, W states and the locally maximally entanglable states. We furthermore show how basic properties of M-spaces can be understood transparently by manipulating their monomial stabilizer groups. In particular, we derive a unified procedure to construct an eigenbasis of any M-space, yielding an explicit formula for each of the eigenstates. We also discuss the computational complexity of M-spaces and show that basic problems, such as estimating local expectation values, are NP-hard. Finally, we prove that a large subclass of M-spaces—containing, in particular, most of the aforementioned examples—can be simulated efficiently classically with a unified method.

and dynamics of stabilizer states in a variety of settings-in fact, the PSF is commonly used in virtually all subfields of quantum information. Important applications include quantum error correction [1], measurement-based computing [2] and classical simulations of quantum circuits [3]. In addition, the PSF is used in condensed matter physics; cf the study of topological order [4].
Notwithstanding its success, a drawback of the PSF is that it describes a relatively small class of states. In particular, there are only finitely many stabilizer states for each given system size. Furthermore, these states have very particular properties: for example, (modulo some trivial cases) they cannot be unique ground states of two-body Hamiltonians, every qubit is either maximally entangled with the rest of the system or completely disentangled from it, most interesting Pauli stabilizer states have zero correlation length and so on [5,6]. In addition, by definition the PSF regards only commuting stabilizers operators. This situation prompts the question of whether it is possible to enlarge the class of stabilizer states while maintaining a transparent stabilizer-type description. Such generalizations could lead to new insights into many-body quantum states as well as novel applications, such as better error-correcting codes, new information-theoretic protocols and new quantum states/processes that can be simulated efficiently classically.
A central feature of the PSF is that relevant information about stabilizer states can be extracted transparently and efficiently by suitably manipulating their stabilizer groups. Aiming at generalizing the PSF, our goal is to identify a mathematical structure that is richer than the PSF while maintaining similarly clearcut maps from the 'stabilizer picture' to the 'state picture'. Our route towards this end starts with the following observation: all Pauli operators are monomial matrices, i.e. precisely one matrix entry per row and per column is nonzero. The basic premise of this work is then to consider arbitrary monomial unitary operators (with efficiently computable matrix elements) as stabilizer operators, giving rise to a general 'monomial stabilizer formalism' (MSF).
Part of the initial motivation for considering the MSF stems from the fact that monomial matrices are simple operators with favorable mathematical properties. In particular, any group generated by a set of monomial matrices consists entirely of monomial matrices, making such groups rather manageable objects. Most of the motivation for studying the MSF, however, is an a posteriori one; in fact one of the purposes of this paper is to argue that the MSF has several nice features and to make a case for the further investigation of this framework as an interesting generalization of the PSF. We will do so by means of the following two contributions (cf section 2 for a more detailed summary).
(a) We show that-perhaps surprisingly-a variety of important quantum many-body states are covered by the MSF, demonstrating that the latter is significantly richer than the PSF. (b) We show that the basic properties of M-spaces can be transparently described by manipulating their monomial stabilizer groups; this will lead, in particular, to efficient classical simulations of a subclass of M-states.
In view of these features and taking into account the wide applicability of the PSF, we believe that the MSF provides a promising avenue to describe and simulate interesting many-body states, potentially leading to new applications. Finally, we refer the reader to an upcoming work [7], where the methods developed in the present work are used to arrive at new, efficient classical simulations of quantum Fourier transforms.

Summary of the results
Next we summarize the main results of this paper, outlining in particular the content of contributions (a) and (b) mentioned above. The +1 common eigenspace of a set of operators will refer to the space of all states which are eigenvectors with eigenvalue 1 of all operators in the set. Generalizing the notion of a Pauli stabilizer code, the +1 common eigenspace of a set of monomial unitary operators will be called an M-space. If an M-space is one-dimensional (1D), its (up to a global phase) unique element will be called an M-state, generalizing the notion of Pauli stabilizer states.
As for (a), we will demonstrate that the MSF encompasses, in addition to all Pauli stabilizer states and codes, the following important state families: the ground level of the Affleck-Kennedy-Lieb-Tasaki (AKLT) model [8]; the ground levels of Kitaev's quantum double models, which describe both Abelian and non-Abelian anyonic systems [4]; the Laughlin wavefunction at filling fraction ν = 1 [9]; the family of locally maximally entanglable (LME) states [10]; coset states of Abelian groups (cf the Abelian hidden subgroup problem [11]); coherent probabilistic computations [12]; W states [13]; and Dicke states [14].
These examples demonstrate the richness of the MSF. They also show that M-states generally do not display the aforementioned 'special' features of Pauli stabilizer states. For example, there are interesting M-states that have non-commuting stabilizer groups and are unique ground states of two-body Hamiltonians.
As for (b), we will establish basic maps from the stabilizer picture to the state picture by showing how a designated orthonormal basis (called here the orbit basis) of any M-space can be constructed when the latter is described in terms of a set of monomial stabilizers. The procedure yields an explicit formula for each basis state, formulated entirely in terms of manipulations on the stabilizer group. This result applies to arbitrary M-spaces and thus, in particular, to all the examples given above. In other words, one obtains a single unified method to treat a number of seemingly unrelated state families. It will also follow from our analysis that all M-states have a common, particularly simple structure, namely the nonzero-amplitudes of any M-state are all equal in modulus.
We will subsequently use the orbit basis construction to investigate classical simulations. Whereas within the PSF many quantities of interest can be computed efficiently for all Pauli stabilizer states, the situation is different for general M-spaces. We will show that one cannot hope for general efficient algorithms for several basic problems (such as estimating local density operators), as we will prove their NP-hardness. In other words, the MSF is too rich a framework to allow for generally applicable efficient simulations. However, it is important that these results regard worst-case complexity. In fact, based on our characterization of the orbit basis, we will identify a relevant subclass of M-states for which efficient classical simulations can nonetheless be achieved; this subclass contains, in particular, almost all the examples listed in (a).
We will illustrate our general constructions with examples throughout this paper. It is noteworthy that our methods allow us to recover, with one unified method, the classical simulatability of a variety of state families, including the Pauli stabilizer states [3] and the quantum double models [15,16]. In addition, the MSF allows us to rederive, in a new and unified way, the standard basis expansion of stabilizer states in terms of cosets of Z 2linear spaces [17] as well as the matrix product state basis of the ground level of the AKLT model [8].
Notations and conventions. All Hilbert spaces considered in this paper are finite dimensional. We will often consider unnormalized quantum states in order not to overload the notation. The group G generated by a set of unitary operators U 1 , . . . , U m (i.e the set obtained by taking arbitrary products of the U i and their inverses) is denoted by G = U 1 , . . . , U m . We also remark that some proofs will be presented in the appendices.

M -states and M -spaces
Let H be a Hilbert space with orthonormal basis where I denotes some finite set. We will mostly consider H to be a multi-party tensor product of d-dimensional local spaces and B will usually be a product basis, but our arguments hold for arbitrary spaces and bases. A unitary operator is called B-monomial if it can be written as a product U = P D, where D is a diagonal unitary operator in the basis B (thus containing phases on its diagonal) and where P is a permutation matrix in this basis. Equivalently, the matrix representation of U in the basis B contains precisely one nonzero entry per row and per column.
If U and U are B-monomial operators, so are the operators U † and UU . Furthermore, if U is B-monomial and V is B -monomial, then U ⊗ V is monomial relative to the tensor product basis B ⊗ B . In the following, it will usually be clear from the context which basis B is considered. Therefore, the prefix 'B-' will mostly be omitted, i.e. we will simply refer to 'monomial' operations. Elements of B will typically be denoted by |x , |y and |z .
An important feature of monomial operations is that products U 1 U 2 · · · U t remain monomial regardless of the length of the product, as long as every U i is monomial (relative to the same basis). This implies, in particular, that if a group is generated by a set of monomial operators, then all elements in this group are monomial.
The +1 common eigenspace of a set {U 1 , . . . , U m } of operators is the space of all |ψ satisfying U i |ψ = |ψ for every i = 1, . . . , m. In this paper, we will only consider finite stabilizer groups G. 2 Otherwise, in the above definitions no restrictions are placed on the operators U i except their unitarity and monomiality. In the following, we envisage M-spaces M which are defined in terms of a set of equations (2). A basic problem will then be to understand how features of M can be traced back to features of G. More stringently, we will study the computational hardness of determining certain quantities (such as expectation values of local observables) by suitably manipulating the generators U i . To make meaningful statements about computational efficiency issues, it needs to be made clear which classical descriptions of the U i are considered to be available (in particular, 6 one should restrict to M-spaces which have efficient descriptions) and how computational cost is measured. This discussion is postponed to section 7. In sections 4-6 we will not yet worry about efficiency issues and our treatment will hold for M-spaces in general.

Examples
The paradigm of M-states and -spaces encompasses a number of important quantum many-body states.
• Kitaev's quantum double models. These are generalizations of the toric code which describe systems of Abelian and non-Abelian anyons [4].
• LME states [10]. This class was recently introduced in studies of multipartite entanglement.
• Coset states of finite Abelian groups. These states occur in important quantum algorithms, namely the Abelian hidden subgroup problem (see, e.g., [11]).
Note that the above examples occur in different areas of application, ranging from topologically ordered systems (cf quantum double models) to quantum algorithms (cf coset states) to multipartite entanglement studies (cf W states and LME states). This richness of M-spaces motivates their study as a general framework, as initiated in this paper. We briefly discuss the Pauli stabilizer states/codes, the AKLT model and the LME states here. See appendix A for a discussion of the other examples.

Pauli stabilizer states and codes
An n-qubit Pauli operator has the form σ = γ σ 1 ⊗ · · · ⊗ σ n , where each σ k is either the singlequbit identity operator or one of the Pauli matrices and where γ = ±1, ±i. A linear subspace of an n-qubit system is a Pauli stabilizer code if there exists a set of commuting Pauli operators that has this space as its +1 common eigenspace. If a Pauli stabilizer code is 1D, then its unique (up to a global phase) element is called a Pauli stabilizer state. Note that the Pauli matrices and the identity are obviously monomial (relative to the standard basis). Since tensor products of monomial matrices are again monomial, we find that every Pauli operator is unitary and monomial relative to the computational basis. Thus every Pauli stabilizer code/state is an M-space/state. Note that general monomial unitary stabilizer groups are significantly richer than the Pauli stabilizer groups. For example, Pauli operators are product operators, whereas elements of general monomial stabilizer groups need not be. Secondly, Pauli stabilizer groups are Abelian, whereas general monomial stabilizer groups can be non-Abelian (cf the AKLT model in 7 section 4.2). Finally, the number of Pauli stabilizer groups (and thus also the number of Pauli stabilizer states/codes) is finite for a given Hilbert space dimension. Monomial unitary stabilizer groups, on the other hand, form a continuous family.
Owing to the richer structure of the allowed stabilizer groups, M-states/spaces may exhibit features that can never be present in Pauli stabilizer states/codes. For example, for every Pauli stabilizer state the bipartite entanglement between any qubit and the rest of the system is either maximal or zero [6]. Also the localizable entanglement [21] between any two qubits is either zero or maximal [6]. It is, however, easy to give examples of M-states where these entanglement measures vary continuously 3 . It is also known that, except for some uninteresting cases, no Pauli stabilizer state can be the unique ground state of a two-body Hamiltonian [5]. The example of the AKLT model (cf section 4.2) shows, however, that interesting M-states with this property do exist.

The Affleck-Kennedy-Lieb-Tasaki (AKLT) model
We show that the ground level of the 1D spin-1 AKLT model is an M-space; to our knowledge this is a new way of describing the AKLT model.
Consider two spin-1 particles with local basis {|0 , |1 , |2 }. Let π be the projector onto the subspace spanned by Now consider a 1D chain of spin-1 particles labeled from 1 to n where one may consider either periodic or open boundary conditions. Let the local two-particle Hamiltonian H i,i+1 act as I − π on spins i and i + 1 and note that H i,i+1 0. Let M denote the zero-energy ground state subspace of the AKLT Hamiltonian H aklt = H i,i+1 [8]. Since H i,i+1 0, the space M coincides with the 0 common eigenspace of the operators H i,i+1 . In the case of open boundary conditions, M is fourfold degenerate; in the case of periodic boundary conditions, M is 1D, i.e. the AKLT ground state is unique [8].
Let U be the following monomial unitary operator: It is straightforward to show that the +1 eigenspaces of U and π coincide. Letting U i,i+1 act as U on spins i and i + 1 it follows that M is the +1 common eigenspace of the U i,i+1 and thus an M-space.

Locally maximally entanglable states
An n-qubit state |ψ is LME if there exists a unitary product operator U = U 1 ⊗ · · · ⊗ U n such that 8 for some complex phases γ x , where |x denote n-qubit computational basis vectors and where the sum is over the entire basis. LME states were recently introduced in [10] where it was shown that |ψ is LME if and only if this state can be maximally entangled to an ancillary n-qubit system using local controlled-unitary operations 4 . It was also shown that the family of LME states contains all Pauli stabilizer states. In turn, all LME states are M-states. This readily follows from the stabilizer description of LMEs introduced in [10], which we briefly repeat here. Define the diagonal unitary operator D := γ x |x x| and let X i be the Pauli X operator acting on qubit i. Since |+ := |0 + |1 is the unique +1 common eigenvector of the n operators X i , it follows that |ψ := D|+ n = γ x |x is the unique +1 common eigenvector of the operators U i := D X i D † . Note that the U i are unitary and monomial relative to the computational basis so that |ψ is an M-state. As |ψ = U |ψ this immediately implies that |ψ is an M-state relative to the product basis B := {U |x }.

Preliminary concepts
Here we introduce some concepts that will be needed in the statement and proofs of our main results below.

The projector ρ
Consider a finite unitary group G with +1 common eigenspace M and let ρ denote the orthogonal projector onto M (remark that here the symbol ρ does not represent a density matrix). Then ρ coincides with the 'group averaging operator' Furthermore, one has For completeness, proofs of these properties are given in appendix B.

The support of a subspace
Consider a Hilbert space H with orthonormal basis B. The B-support, or simply the support, of a state |ψ is the set of basis states |x ∈ B for which x|ψ is nonzero. The support of a linear subspace M of H is the union of the supports of its elements. The support of a state/subspace will be denoted by supp(|ψ)) and supp(M), respectively. Letting ρ denote the orthogonal projector on the space M, it is easily verified that the following statements are equivalent.

Uniform superpositions
Consider a subset O ⊆ B and a complex phase γ y for every |y ∈ O. Consider the states The state |ψ is called a B-uniform superposition, or simply a uniform superposition. Note that O is precisely the support of |ψ . The state |O is called the equal superposition over O.

The permutation group P
Consider a unitary monomial operator U = P D where P is a permutation and D is diagonal; note that this decomposition is unique. Henceforth we will denoteŪ := P. Let G = U 1 , . . . , U m be a unitary monomial group. The set P = Ū : U ∈ G is called the permutation group of G. Using that U V =ŪV for every unitary monomial U and V , it can easily be shown that P is a group generated by the operatorsŪ i . 5 The group P naturally acts as a permutation group on B. Consider the orbit of |x : It is well known that every |x belongs to precisely one orbit (i.e. the orbit O x ). Straightforwardly applying the definition of P one finds that |y ∈ O x if and only if there exists a complex phase ξ and U ∈ G such that U |x = ξ |y . This elementary property will often be used in the following.
Example. Throughout this section as well as section 6, we illustrate the various concepts with a simple 'running example'. Consider the n-qubit diagonal operator Further, let S i denote the SWAP gate acting on qubits i and i + 1, where i = 1, . . . , n − 1. The operators T and S i are unitary and monomial (relative to the computational basis). We consider the group G w generated by these operators. Since T is diagonal, one hasT = I . Since S i is a permutation matrix, one hasS i = S i . Thus P w is generated by the operators S i . It is easily verified that P w has n + 1 orbits O 0 , . . . , O n , where O i contains all computational basis states |x , where x is a bit string with Hamming weight i (that is, precisely i entries are equal to 1).

The phases ξ x (y)
Let M be an M-space with stabilizer group G and let ρ be the orthogonal projector on M. Fix |x ∈ supp(M) (hence ρ|x is nonzero) and |y ∈ O x arbitrarily. Then there exists U ∈ G and a complex phase ξ such that U |x = ξ |y . We consider the set x→y of all phases that may occur in this way: Remarkably, this set is in fact a singleton. To see this, consider ξ ∈ x→y and a corresponding U ∈ G. Using that ρ = ρU it follows that ρ|x = ξρ|y . Thus ρ|y is proportional to ρ|x and the proportionality factor is precisely given by ξ . This shows that this phase does not depend on U , so that x→y is indeed a singleton. In the following, we will denote the unique element of this set by ξ x (y). We have proved: For every |x ∈ supp(M) and |y ∈ O x , one has ρ|x = ξ x (y)ρ|y .
We emphasize that the phase ξ x (y) is only defined for |x ∈ supp( M) and |y ∈ O x . Finally, note that the phase ξ x (x) is well defined for every |x in the support of M, since |x ∈ O x . Using that I ∈ G and I |x = |x , we in fact find that ξ x (x) = 1.
Example. Consider the group G w as above and let M w be its +1 common eigenspace. We show in section 6 that the basis state |e 1 := |10 · 0 belongs to the support of M w . Note that |e 1 belongs to the orbit O 1 . The other vectors in this orbit are |e 2 , . . . , |e n where e i denotes an n-bit string with a 1 in the ith slot and zeros elsewhere. Letting P i denote the operator which swaps qubits 1 and i (which can easily be obtained as a suitable product of S k gates), one has P i |e 1 = |e i . It follows that ξ e 1 (e i ) = 1 for every i.

The orbit basis
Consider an M-space M specified in terms of a stabilizer group G. The latter may, for example, be given in terms of a set of generators. Our goal is to construct an orthonormal basis of M assuming no prior information about this space except for the group G. In this section, we prove two results which will directly lead to such a construction. The first of these theorems characterizes the support of M, the second theorem characterizes a designated basis of M called the orbit basis.

Characterizing the support
We need some notation. For every |x ∈ B, let G x be the set of all U ∈ G which have |x as an eigenvector. This set is a subgroup of G. We let {U x,1 , . . . , U x,l } be an arbitrary set of generators 6 of G x .

Theorem 1 (support of M-space). Consider an M-space
Then the following statements are equivalent.
Proof. Lemma 1 shows that for every |x in the support of M and |y ∈ O x , one has ρ|y = 0 so that also |y ∈ supp(M). This shows that the entire orbit O x must be contained in the support. Thus supp(M) is a union of orbits as claimed.
[a⇒c] If O x ⊆ supp(M), then |x belongs to the support of M since |x ∈ O x . For every |x in the support, the phase ξ x (x) is well defined and in fact ξ x (x) = 1 (cf section 5.5).
Consider an arbitrary U ∈ G x , i.e. there exists a phase ξ such that U |x = ξ |x . By uniqueness of the phase ξ x (x) it follows that ξ = ξ x (x) = 1.
[c⇒b] Consider an arbitrary U ∈ G. Then U |x = ξ |y for some complex phase ξ and some |y . If y = x, then x|U |x is zero. If y = x, then U ∈ G x and thus ξ is equal to +1 by assumption (c). Thus x|U |x = ξ = 1.
[b⇒a] Let ρ denote the orthogonal projector onto M. Using (6) we have Owing to (b), every term x|U |x in the sum is non-negative. Moreover, there is at least one nonzero term, i.e. when U is the identity. This shows that x|ρ|x is nonzero so that |x belongs to the support of M. But then O x ⊆ supp(M) since the support is a union of orbits.
[c⇒d] This equivalence is straightforward since the U x,i generate G x .
Example. Consider the M-space M w with stabilizer group G w as in section 5. For every computational basis state |x , one has where |x| denotes the Hamming weight of x. This shows that, for every x with |x| = 1, one has T |x = λ|x for some λ = 1. Invoking theorem 1, we find that none of the orbits O i with i = 1 are contained in the support of M w . On the other hand, the orbit O 1 is contained in the support.
To show this, consider the n-qubit W state It is straightforward to verify that T |W = |W = S i |W for every i, showing that |W ∈ M w . Since this state has support O 1 , it follows that O 1 ⊆ supp(M) w . In conclusion, the support of M w is identical to O 1 .

Constructing the orbit basis
For every |x ∈ supp(M), consider its normalized projection onto M: Henceforth whenever considering a state (15) we will tacitly assume that |x belongs to the support since otherwise ρ|x = 0. Note that by construction U |ψ x = |ψ x for every U ∈ G. Interestingly, the states |ψ x have the following explicit form.

Lemma 2 (orbit states).
Every |ψ x is a uniform superposition state given by Proof. Using that ρ is given by (6), one has 12 for some coefficients c y . This shows that supp(|ψ x ) ⊆ O x . Furthermore, for every |y ∈ O x there exists U ∈ G such that U |x = ξ x (y)|y . Using that U † |ψ x = |ψ x , it follows that In combination with (17) this implies that Since |ψ x is normalized, it follows that: for some complex phase γ . Finally, it follows from definition (15) that x|ψ x = ρ|x > 0. This shows that γ = 1.
Since |ψ x is a uniform superposition over an orbit of P, we call this state an orbit state. More precisely we call |ψ x the orbit state determined by |x .
Example. The orbit state of M w determined by |e 1 is the n-qubit W -state: Next we show that a basis of M can be constructed by selecting a suitable subset of orbit states; this is the orbit basis. Proof. Since B is a basis of the Hilbert space and since |ψ x is defined as the projection of |x onto M, the collection of all orbit states spans M (although generally these states are not linearly independent). Consider an arbitrary |x in the support of M. By construction there exists an i between 1 and d such that |x ∈ O i . Lemma 1 shows that ρ|x ∝ ρ|x i so that |ψ x ∝ |ψ i . Since the orbit states span M, it follows that the states |ψ i span M as well. Since |ψ i has the orbit O i as its support owing to (16) and since these orbits are mutually disjoint, the states |ψ i are mutually orthogonal. This shows that is an orthonormal basis. Finally, lemma 1 shows that the basis is independent (up to global phases) of the choice of orbit representatives |x i . Theorems 1 and 2 show that the following procedure allows us to correctly identify the orbit basis of M.
• Determine all orbits of P and consider a representative x k in each orbit O k .
• For each k, decide whether O k ⊆ supp(M) by means of the characterization in theorem 1.
• The orbit basis is the collection of all orbit states |ψ x k for which O k ⊆ supp(M).
As desired, this procedure can be implemented by means of manipulations on the stabilizer group G. Furthermore, each orbit state is itself characterized completely in terms of properties of G owing to lemma 2.
Example. We have seen that supp(M w ) = O 1 . Theorem 2 thus shows that M w is 1D, with |ψ e 1 as its unique element. Furthermore, we have shown that |ψ e 1 = |W . Thus, the W state is an M-state with stabilizer group G w .
Note that the orbit basis construction applies to arbitrary M-spaces and thus, in particular, to all instances given in section 4. This yields one unified method to analyze all these state families; see section 8 for further illustrations.

Some corollaries
Next we discuss some immediate corollaries of theorems 1 and 2. The first corollary is interesting in that it bounds the dimension of M by means of a purely combinatorial quantity, namely the number of orbits. Note that there exist well-developed tools to count/estimate orbits of permutation groups, which may thus be imported into the study of M-spaces.

Corollary 1 (dimension). The number of orbits d contained in the support of M (cf theorem 1) coincides with the dimension of M. It follows that this dimension is upper bounded by the total number of orbits of P.
The above corollary will be used in section 8 to give a simple proof that the AKLT ground level in the open boundary conditions case is fourfold degenerate.
Secondly, theorem 2 leads to the following characterization of M-states.

Corollary 2 (M-states).
Every M-state |ψ is a uniform superposition. More precisely, if |x is an arbitrary basis state satisfying x|ψ = 0, then |ψ ∝ |ψ x where |ψ x is given explicitly in lemma 2.
Note that the support of any M-state coincides with precisely one orbit. This implies in particular that knowledge of a single |x satisfying x|ψ = 0 implies complete knowledge of the entire support of |ψ , which then must coincide with O x . Theorem 1 may be used to determine a suitable x such that |ψ ∝ |ψ x .
It is interesting to compare corollary 2 with a characterization of Pauli stabilizer states obtained in [17]. Consider a Pauli stabilizer state |ψ on n qubits. Then there exists a linear subspace S of Z n 2 , an x ∈ Z n 2 and complex phases ξ y such that Thus every Pauli stabilizer states is a uniform superposition, the support of which is identified with a coset x + S. Corollary 2 shows that, in fact, all M-states have a basis expansion with an analogous structure. Corollary 2 can in fact be used to rederive (22) in a simple way; see section 8.
Restricting attention to many-qubit systems there is a noteworthy connection to LME states. In the terminology of this paper, an n-qubit state |ψ is LME iff there exists a product basis B such that this state is a B-uniform superposition where supp(|ψ ) is the entire basis B. Corollary 2 shows that n-qubit M-states (when monomiality is considered relative to product bases) can be regarded as generalizations of LME states where uniform superpositions with arbitrary supports are considered.
Finally, we discuss a subclass of stabilizer groups for which the results above can be simplified. A monomial unitary group G is called pure if it has a generating set of the form {P 1 , . . . , P k , 1 , . . . , m } where every P i is a permutation and every j is diagonal. Such a generating set is also called pure. For every |x , let |O x denote the equal superposition over the orbit O x (recall section 5.3).

Corollary 3 (pure stabilizer groups). Let M be an M-space with pure stabilizer group G.
Then the orbit basis of M has the form = {|O x 1 , . . . , |O x d }.
Proof. Let P be the permutation group of G. Consider a pure generating set {P i , j } of G. Sincē P i = P i and¯ j = I , it follows that P is generated by the operators P i . This implies that P ⊆ G.
Let |x belong to the support of M and consider an arbitrary |y ∈ O x . Then there exists P ∈ P such that P|x = |y . Since P ∈ G and owing to the uniqueness of the coefficient ξ x (y) it follows that ξ x (y) = 1. Invoking lemma 2 it follows that |ψ x = |O x .
Examples of pure stabilizer groups are numerous: consider, e.g., the toric code states, quantum double models, W -states, coherent probabilistic computations and coset states of Abelian groups (cf section 4 and appendix A). We will consider the example of quantum doubles in more detail in section 8.

Computational complexity and classical simulations
So far, we have studied general mathematical features of M-spaces, not worrying about the computational complexity of the concepts involved. Here we discuss such issues. For concreteness, we consider many-qubit systems where monomiality is defined relative to the computational basis {|x }; generalizations are straightforward.

Classical descriptions of M-spaces
First, it needs to be made clear which classical descriptions of M-spaces are considered to be available. It is natural to consider n-qubit M-spaces M that meet the following requirements.
(i) M has a stabilizer group G = U 1 , . . . , U m with m = poly(n) generators. A classical description of each U i is considered to be given as an input; each of these descriptions is assumed to be efficient. (ii) Every generator U i has efficiently computable matrix elements in the following sense.
Suppose that U i acts on the computational basis as for some complex phases λ i (x) and for some permutation π i of the set of n-bit strings. Then U i is said to be efficiently computable if, given any n-bit string x as an input, the following two conditions are fulfilled 7 .
-There exists a poly(n) time classical algorithm to compute π i (x) and π −1 i (x). -There exists a poly(n, k) time classical algorithm to compute the phase λ i (x) up to k bits of precision. These conditions entail that, given any row of the matrix U i , it is possible to efficiently compute which matrix element within that row is nonzero and what the value of that matrix element is, and a similar condition for the columns.
The above conditions are met in many cases of interest. In particular, we have:

except for the LME states, all M-spaces considered in section 4 satisfy (i)-(ii).
This statement is easily verified and the arguments are omitted here. LME states generally do not satisfy (i)-(ii) since the phases γ x in (5) may not be efficiently computable. However, all LME states where the function x → γ x is classically computable in poly(n, k) time up to k bits-which is a large class-do satisfy conditions (i)-(ii).

Computational complexity
Considering n-qubit M-spaces described as above one may investigate the computational hardness of a variety of tasks; here an algorithm is considered to be efficient if it runs in poly(n) time. Natural problems are as follows. We will show that one cannot hope for efficient classical algorithms for P1-P3 that apply to all M-spaces in general. This is a point where the MSF sharply contrasts with the PSM, where many problems of interest can be answered efficiently for arbitrary Pauli stabilizer states and codes; this holds in the particular problems P1-P3. Hardness of the above problems shows that one should look for relevant subclasses of M-spaces for which efficient solutions are possible. The intractability of P1-P3 in fact holds even for very simple M-spaces, namely those with diagonal, local generators U i . Problem 1. The input is a set of n-qubit diagonal unitary operators {U 1 , . . . , U m }. Each U i acts nontrivially on at most three qubits and all matrix entries of U i (in the computational basis) are 0, 1 or −1. The problem is to decide whether these operators have a +1 common eigenvector.

Problem 2.
The input is a set of diagonal unitary operators U i as in problem 1, with the following additional constraint: it is promised that these operators have a unique (up to a global phase) +1 common eigenstate |ψ , which is thus an M-state. The problem is to sample classically from the distribution {| x|ψ | 2 }.

Problem 3.
The input is as in problem 2. The problem is to compute ψ|Z i |ψ with accuracy = 1/poly(n), where Z i is the Pauli σ z operator acting on qubit i.

Theorem 3.
None of the problems 1-3 can be solved classically in polynomial time unless P = NP.
Theorem 3 will be proved by reductions to (variants of) the satisfiability problem, which is NP-complete. See appendix C.

Classical simulations
We now focus on the classical simulation problems P2 and P3 in more detail. Even though these tasks are intractable in their worst case, efficient solutions exist for subclasses of M-states. Here we provide sufficient criteria for the existence of efficient algorithms. Then there exists an efficient classical algorithm to sample the distribution := {| y|ψ | 2 }. Furthermore, then there exists an efficient classical algorithm to estimate ψ|A|ψ with accuracy = 1/poly(n) with success probability that is exponentially close to 1.
Proof. Since |ψ is a uniform superposition with support O x for some |x (recall corollary 2), is the uniform distribution over this orbit. This shows that efficient classical algorithms for (a) and (b) imply an efficient algorithm to sample . We consider the second claim. Every klocal observable can be written as a linear combination of poly(n) Pauli operators where each coefficient in the linear combination has modulus not greater than 1. Therefore it suffices to prove the claim for Pauli operators. For every a ∈ Z n 2 , consider X (a) := X a 1 ⊗ · · · ⊗ X a n , Z (a) := Z a 1 ⊗ · · · ⊗ Z a n .

(24)
Then every Pauli operator can be written as γ X (a)Z (b) for some a, b and some γ ∈ {±1, ±i}. See appendix E for some standard properties of Pauli operators. For every |y ∈ O x , define Using that |ψ ∝ |ψ x and lemma 2, one finds that By a standard Chernoff bound argument (see, e.g., the appendix of [23]), the sum in (26) may be estimated with 1/poly(n) error with exponentially small failure probability by generating poly(n) random elements |y i ∈ O x and by computing the average of F(y i ). Furthermore, owing to assumptions (a)-(d) this procedure can be implemented in polynomial time.
Note that (a) can be approached with theorem 1. Furthermore, in order to compute ξ x (y), it suffices to determine any single U ∈ G such that U |x ∝ |y , since then ξ x (y) is simply given by the matrix element y|U |x (cf section 5.5). Finally, we point out the following elementary approach to (b).

Lemma 3.
The following procedure generates a random |y ∈ O x . First, generate a random permutation P ∈ P. Then compute |y = P|x and output y.
The proof of lemma 3 uses basic group theory arguments and is given in appendix D. Interestingly, there exists a well-developed theory of (approximately) generating random elements in finite groups and efficient algorithms are available for a variety of groups (see, e.g., [24]). These methods may be thus imported into the study of classical simulations of M-states.
Efficient algorithms for (a)-(d) exist for a variety of M-states. In fact it can be shown that:

except for the LME states, all M-states considered in section 4 satisfy (a)-(d) in theorem 4.
As before, LME states for which the function x → γ x can be computed efficiently do satisfy (a)-(d). In section 8, we work out the examples of stabilizer states and quantum double models.

Applications
Throughout this paper we have illustrated our results by means of the simple example of the W states. Here we give some more sophisticated examples, namely the Pauli stabilizer states, the AKLT model and Kitaev's quantum double models.

Pauli stabilizer states
We show how the MSF can be used to rederive some interesting features of Pauli stabilizer states in a new way. In particular: (i) we rederive the expansion (22) first proved in [17]; (ii) we show that the conditions (a)-(d) of theorem 4 are fulfilled, thus showing that Pauli stabilizer states can be efficiently simulated classically in the sense of theorem 4; this recovers (a variant of) the result [3].
Let |ψ be an n-qubit Pauli stabilizer state with Pauli stabilizer group G. Recall that every minimal set of generators of G contains precisely n Pauli operators [22]. We arbitrarily fix such generators {σ 1 , . . . , σ n }. Since Y = i X Z , every generator can be written as where γ i ∈ {±1, ±i}; recall also the notation (24). Remark that σ i is a product of a permutation matrix X (s i ) and a diagonal matrix γ i Z (t i ) so thatσ i = X (s i ). Let S be the Z 2 -linear space generated by the vectors s i . Using (E.1) and (E.3) one easily finds that for every x ∈ Z n 2 . Invoking corollary 2 these identities immediately imply that |ψ has the form (22) for some x ∈ Z n 2 . We now address (ii). First, we show that an x such that |ψ ∝ |ψ x can be computed efficiently. Recall the definition of the group G x given in section 6.1. Let D denote the subgroup consisting of all σ ∈ G satisfying σ ∝ Z (b) for some b. It follows from (E.1) and (E.2) that G x = D for every x. We now claim: generating set {D 1 , . . . , D l } of D can be determined in poly-time; moreover, D j = (−1) u j Z (d j ) for some (efficiently computable) u j ∈ Z 2 and d j ∈ Z n 2 . This result can be proved using standard Pauli stabilizer arguments; for completeness a proof is given in appendix E. Theorem 1 now implies that O x = supp(|ψ ) if and only if D i |x = |x for every i. Using (E.2) this is equivalent to requiring that x T d j = u j for every j. A solution x to this system of equations can be computed efficiently.
Since we have access to a generating set {s 1 , . . . , s n } of S, we can efficiently determine whether y ∈ S + x, given y as the input. Also a random element in x + S can easily be generated efficiently.
Finally, we show that, given any y ∈ x + S, the phase ξ x (y) can be computed efficiently. First we compute an arbitrary a = (a 1 , . . . , a n ) ∈ Z n 2 satisfying a i s i = x + y; this regards solving a system of linear equations and can be done efficiently. Using the properties (E.1)-(E.5) it follows that for such a one has It follows that σ a 1 1 · · · σ a n n |x = ξ |y for some complex phase ξ ; by the uniqueness of ξ x (y) we have ξ = ξ x (y). Given a, x and y it is straightforward to compute ξ in poly-time.

AKLT model with open boundary conditions
Consider the AKLT model with open boundary conditions. Let M open denote the ground level subspace. We will use the MSF to prove the following properties (see also [8,25]). (i) The ground level is fourfold degenerate. (ii) An orthonormal basis of ground states is given by |ψ σ open = Tr{σ σ a 1 · · · σ a n }|a 1 · · · a n .
We prove the results for n even. Let G open be the group generated by the operators U i,i+1 . Let P denote the permutation operator obtained by replacing all minus signs in (4) by plus signs. Then obviouslyŪ i,i+1 = P i,i+1 so that P open is generated by the operators P i,i+1 . If n is even, then it is straightforward to show that P open has four orbits.
• O I contains all basis states with an even number of 0s, an even number of 1s and an even number of 2s. A representative element is |a I = |0 · · · 000 .
• O X contains all basis states with an even number of 0s, an odd number of 1s and an odd number of 2s. A representative element is |a X = |0 · · · 012 .
• O Y contains all basis states with an even number of 1s, an odd number of 0s and an odd number of 2s. A representative element is |a Y = |1 · · · 102 .
• O Z contains all basis states with an even number of 2s, an odd number of 0s and an odd number of 1s. A representative element is |a Z = |2 · · · 201 .
The Pauli matrices satisfy the commutation relations By using the commutation relations (30) and by applying the definition of the phases ξ x (y), it can also be shown directly that |ψ σ open is the orbit state determined by |a σ . This argument is omitted here.

AKLT model with periodic boundary conditions
Consider the AKLT model with periodic boundary conditions and n even. Let M per denote the ground level subspace. We rederive the following properties, first proved in [8].
Using the commutation relations (30) it is straightforward to show that |ψ per belongs to M per . As above, P per is generated by the permutation matrices P i,i+1 (now imposing periodic boundary conditions). Furthermore, one may verify that P per has the same four orbits as in the open boundary conditions case. Note also that |ψ per has support O I so that O I ⊆ supp(M per ). We show that none of the other orbits are contained in the support. To do so define the operators A = U n,n−1 U n−1,n−2 · · · U 2,1 U 1,n , One can then verify that B|a σ = −|a σ for every σ = X, Y, Z . Since B ∈ G per , owing to theorem 1 this implies that none of the orbits O X , O Y or O Z are contained in the support of M per . It follows that the support coincides with the single orbit O I . Invoking corollary 1 then implies that M per is 1D.

Quantum double models
Consider a quantum double model defined on a sphere (cf appendix A.2). Such a system has a unique ground state |ψ qd [4], which is thus an M-state. Letting G be the stabilizer group of this state, here we show: G meets requirements (a)-(d) of theorem 4, showing that |ψ qd can be simulated classically in the sense of theorem 4.
The work [4] can also be used to show that standard basis measurements can be simulated classically (although it is not explicitly discussed there). Efficient simulations of local expectation values were previously achieved using tensor network methods [15,16]. First we describe the stabilizer group. Recall the definition of the operators U p and V v (k) defined in appendix A.2. Since U p is diagonal and V v (k) is a permutation matrix, the group G generated by these operators is pure. Note that the U p mutually commute since these are diagonal operators. Furthermore, it is easily verified that [V v (k), U p ] = 0 and [V v (k), V w (l)] = 0 for all vertices v = w, for every plaquette p and for every k, l ∈ G. It follows that a general element U ∈ G has the form U = P D where Since P is a permutation matrix and D is diagonal, one hasŪ = P.
Next we determine the support of |ψ qd . A standard basis state has the form |g = |g e where the product is over all edges e and where g e ∈ G. Let S be the set of all |g satisfying U p |g = |g for all plaquettes p. We claim that supp(|ψ qd ) = S. To see this, first note that theorem 1 shows that supp(|ψ qd ) ⊆ S. To prove the reverse inclusion, consider an arbitrary U = P D ∈ G as in (33). If |g ∈ S, then D|g = |g . Consequently, where the last inclusion holds since P is a permutation matrix. Invoking theorem 1 proves the claim.
Recall that |ψ qd is an M-state and that G is pure. Corollary 3 implies that |ψ qd is the equal superposition over S and that S must be an orbit of P. Since one manifestly has |e ∈ S (where e has the neutral element e in all its entries), it follows that S = O e . Consequently, |ψ qd coincides with the orbit state |ψ e . This shows that (a) in theorem 4 is fulfilled.
To show (b), we use lemma 3. Since an arbitrary P ∈ P has the form given in (33), a random P can efficiently be generated by generating a random k v ∈ G for every vertex v. Applying P to |e yields a random state in O e .
Condition (c) holds as the support O e = S is defined in terms of polynomially many constraints U p |g = |g , each of which is easily verified. Finally, since |ψ qd is an equal superposition, one has ξ e (g) = 1 for every g ∈ O e .

Outlook
There are several interesting questions to consider for further research. For example, so far we have analyzed existing state families with the machinery of monomial stabilizers; using these techniques it should, however, also be possible to construct and study new interesting M-states and -spaces. For example, can one generalize the present analysis of the 1D AKLT model to arrive at interesting 2D systems? It would also be worthwhile to study suitable parent Hamiltonians for M-states/spaces, the gaps of such Hamiltonians, the relation between the MSF and matrix product states and, more generally, projected-entangled-pair-states, etc. As a different aspect, it would be interesting to obtain a physical interpretation (if it exists) for M-states/spaces (e.g. along the lines of the definition of LME states [10]). Investigating the possibility of efficiently preparing M-states on a quantum computer would be another interesting problem.
where δ(e) = 1 and δ(k) = 0 for every e = k ∈ G (where e is the neutral element of G). To define A v consider the four edges ( f 1 , f 2 , f 3 , f 4 ) incident on v. The signature t i of f i is now +1 if f i is pointing towards v and −1 if it is pointing away from v. Consider a basis state |g 1 , g 2 , g 3 , g 4 where g i ∈ G is associated with the system on edge f i . If k ∈ G define The ground state subspace M qd of H qd consists of those |ψ which are joint +1 eigenvectors of all projectors B p and A v . To show that M qd is an M-space, define where k ∈ G. Then U p is a diagonal unitary operator and V v (k) is a permutation operator. It is immediate that the +1 eigenspaces of U p and B p coincide. Furthermore, the +1 eigenspace of A v coincides with the +1 joint eigenspace of the set To see this, note that (a) A v is a group and (b) A v is the averaging operator over this group. It follows that (cf appendix B) A v is the projector onto the +1 common eigenspace of A v . In other words, the +1 eigenspace of A v coincides with the +1 joint eigenspace of A v .
In conclusion, M qd coincides with the +1 common eigenspace of the operators U p and V v (k) where p ranges over all plaquettes, v ranges over all vertices and k over the entire group G. This shows that M qd is an M-space.

A.3. Laughlin wavefunction at ν = 1
The Laughlin wave function for n particles and filling fraction ν = 1 is |ψ l = n−1 a 1 ,...,a n =0 a 1 ···a n |a 1 · · · |a n . (A.6) Here |a i is a wavefunction for particle i given by where z is a complex number encoding the position of the particle in a 2D plane. Furthermore, is the Levi-Civita completely antisymmetric tensor in n dimensions. Let S i j be the operator which swaps systems i and j. It is straightforward that each operator −S i j is unitary and monomial and that |ψ l is the unique +1 common eigenvector of these operators. Thus |ψ l is an M-state.

A.4. Coherent probabilistic computations
Consider a poly-size classical circuit composed of reversible gates (e.g. Toffoli gates and NOT gates). The circuit acts on an n-bit input string where, say, the first k bits are uniformly random and the last n − k input bits are initialized in 0. Let π x : x ∈ Z n 2 denote the output probability distribution over the set of n-bit strings and define the state |π = √ π x |x . LetC denote the natural translation of C into an n-qubit quantum circuit. Furthermore, let X i and Z j denote the Pauli X and Z operators acting on qubit i and j, respectively, and define P i =C X iC † and D j =C Z jC † . These operators are unitary and monomial (relative to the computational basis) and have |π as unique +1 common eigenstate. Thus |π is an M-state.

A.5. Coset states of finite Abelian groups
Consider a Hilbert space with basis {|g : g ∈ G} where basis vectors are labeled by the elements of a finite Abelian group G. Let χ g be the character of G canonically associated with g and define the operators X (g) : |x → |x + g ; Z (g) : |x → χ g (x)|x , (A.8) for every g, x ∈ G. Note that these operators are unitary and monomial. Consider a subgroup H of G with generating set {h 1 , . . . , h n }. The dual group H ⊥ is the set of all k satisfying χ h (k) = 1 for every h ∈ H ; consider a generating set {k 1 , . . . , k m }. Now define the 'coset state' |H + x = |x + h where the sum is over all h ∈ H . This state is the unique +1 common eigenvector of the operators X (h i ) and Z (k j ) where i ranges from 1 to n and where j ranges from 1 to m. Thus such coset states are M-states. Remark that these states play an important role in the context of quantum algorithms, namely the Abelian hidden subgroup problem.

A.6. W-states and Dicke states
Consider a system of n qubits with computational basis B. Fix k ∈ {1, . . . , n − 1} and define |D k to be the (properly normalized) equal superposition over all basis states |x for which the bit string x has hamming weight k. Note that |D 1 = |W is the W state. Denote T k = α −k · ⊗n , where = diag(1, α) and α = e 2πi n . Further, let S i denote the SWAP gate acting on qubits i and i + 1. The operators T k and S i are unitary and monomial and have |D k as unique +1 common eigenstate. Thus |D k is an M-state.

Appendix B. The projector ρ
Let π = |G| −1 U be the averaging operator over G. We prove that π = ρ. Since π is defined as a sum over all elements in G, one has U π = π for every U ∈ G. Consequently, Moreover, π = π † since for every U ∈ G one has U † ∈ G as G is a group. Thus π 2 = π = π † , i.e. π is an orthogonal projector. Now let W denote the +1 eigenspace of π . First note that the definition of π immediately implies that π |ψ = |ψ for every |ψ ∈ M, showing that M ⊆ W. Furthermore, consider an arbitrary |ϕ ∈ W, i.e. π |ϕ = |ϕ . Since U π = π for every U ∈ G, it follows that: showing that |ϕ ∈ M. This proves that W ⊆ M and we conclude that W = M. Thus ρ and π are orthogonal projectors on the same space, so that these operators are equal.
To prove (7), consider an arbitrary fixed U ∈ G. Since G is a group, the sum V ∈G U V is precisely the sum over all elements of G. Using (6) it follows that Uρ = ρ. If σ = i k X (s)Z (t) (where k ∈ {0, 1, 2, 3} and s, t ∈ Z n 2 ) is a Pauli operator, we will call (k, s, t) the label of σ .

Proof of lemma 4:
We will need to take into account the following technicality. Remark that σ 2 = α I for every Pauli operator σ , where γ = ±1 . Now suppose that there exists σ ∈ G where γ is not equal to 1. Then, since γ I = σ 2 ∈ G one would have γ I |ψ = |ψ so that |ψ = 0, leading a contradiction. Thus we conclude that σ 2 = I for every σ ∈ G. Since the σ i mutually commute and square to the identity, every operator in G can be parameterized as σ (a) := σ a 1 1 · σ a n n , where a = (a 1 , . . . , a n ) ∈ Z n 2 . Using (E.1-E.5) it follows that σ (a) ∝ X ( a i s i )Z ( a i t i ).
(E. 6) This implies that σ (a) ∈ D if and only if a solves the linear equation a i s i = 0. Let {a 1 , . . . , a k } denote a basis of solutions; such a basis can be computed efficiently. Thus σ (a) ∈ D if and only if a = y j a j for some y j ∈ Z 2 . Furthermore, using the definition of σ (a) it is straightforward to show that σ ( y j a j ) = σ (a 1 ) y 1 · · · σ (a k ) y k . (E.7) This shows that the operators D j := σ (a j ) form a generating set of D. Since σ (a j ) is a product of at most n Pauli products, the label of each σ (a j ) can be computed in poly(n) time. Finally, since every element of G squares to the identity, it follows that D j = (−1) u j Z (d j ) for some u j ∈ Z 2 and some d j ∈ Z n 2 .