Two-local qubit Hamiltonians: when are they stoquastic?

We examine the problem of determining if a 2-local Hamiltonian is stoquastic by local basis changes. We analyze this problem for two-qubit Hamiltonians, presenting some basic tools and giving a concrete example where using unitaries beyond Clifford rotations is required in order to decide stoquasticity. We report on simple results for $n$-qubit Hamiltonians with identical 2-local terms on bipartite graphs. Our most significant result is that we give an efficient algorithm to determine whether an arbitrary $n$-qubit XYZ Heisenberg Hamiltonian is stoquastic by local basis changes.


Introduction and Motivation
The notion of stoquastic Hamiltonians was introduced in [1] in the context of quantum complexity theory. The definition of stoquastic Hamiltonians aims to capture Hamiltonians which do not suffer from the sign problem: the definition ensures that the groundstate of the Hamiltonian has nonnegative amplitudes in some local basis while the matrix exp(−H/kT ) is an entrywise nonnegative matrix for any temperature T in this basis. Stoquastic Hamiltonians are quite ubiquitous. Any (mechanical) Hamiltonian defined with conjugate variables [x i ,p j ] = i δ ij of the form with invertible (mass) matrix C > 0, is stoquastic, since the kinetic term is off-diagonal in the position x-basis.
The paradigm of circuit-QED produces such Hamiltonians by the canonical quantization of an electric circuit represented by its Lagrangian density. For the subclass of stoquastic Hamiltonians the problem of determining the lowest eigenvalue might be easier than for general Hamiltonians. In complexity terms the stoquastic lowest-eigenvalue problem is StoqMAcomplete while it is QMA-complete for general Hamiltonians.
This difference in complexity goes hand-in-hand with the existence of heuristic computational methods for determining the lowest eigenvalue of stoquastic Hamiltonians using Monte Carlo techniques [2,3]. These methods can be viewed as stochastic realizations of the power method, namely a stochastic efficient implementation of the repeated application of the nonnegative matrix exp(−H). However, there is no proof that these heuristic methods are as efficient as quantum phase estimation for estimating ground-state energies (some results are obtained in [4]).
Beyond the lowest eigenvalue problem, one can also consider the complexity of evaluating the partition function Z = Tr [exp(−H/kT )] of a stoquastic Hamiltonian. The nonnegativity of exp(−H/kT ) directly ensures that one can rewrite Z as the partition function of a classical model in one dimension higher, which can then be sampled using classical stochastic techniques. Again, rigorous results on this well-known path-integral quantum Monte Carlo method (see e.g. [5]) are sparse: [6] shows how to estimate the partition function with kT = O(log n) for 1D stoquastic Hamiltonians using a rigorous analysis of the path-integral Monte Carlo method.
Stoquastic Hamiltonians play an important role in adiabatic computation [7]. It has been proven that there is no intrinsic quantum speed-up in stoquastic adiabatic computation which uses frustration-free stoquastic Hamiltonians: a classical polynomial-time algorithm for simulating such adiabatic computation was described in [8]. The quantum power of the well-known quantum annealing method which uses the transverse field Ising model is still an intriguing open question, given the road-block cases that have been thrown up for the path-integral Quantum Monte Carlo method in [9] as well as the diffusion Monte Carlo method in [10]. However, given that stoquastic adiabatic computation such as quantum annealing is amenable to heuristic Monte Carlo methods, research has also moved into the direction of finding ways of engineering non-stoquastic Hamiltonians (see e.g. [11], [12], [13]).
If it is true that questions pertaining to stoquastic Hamiltonians are easier to answer than those pertaining to general Hamiltonians, then it is clearly important to be to able to recognize which Hamiltonians are stoquastic.
Work aimed at recognizing where and when a sign problem does not exist, is not new, but few systematic approaches to this problem exist. In this paper we report on our first results in this direction (see e.g. [14]), summarized in Section 1.2.
The definition of being stoquastic is basis-dependent, but it is motivated by computational efficiency. The easiest manner in which one can use the basisdependence is to consider a change of basis implemented by local single-qudit or qubit rotations. This basis change can be efficiently executed in that the description of the k-local Hamiltonian in this new basis can be easily obtained. This paper will focus largely on single qubit rotations. It is however clear that the subclass of single-qubit unitary transformations are only the simplest example of a mapping or encoding of one Hamiltonian H to another Hamiltonian H sim as has been formalized in [15]. Such a more general mapping could allow: constant-depth circuits; ancilla qubits; perturbative gadgets and approximations such that only λ(H) ≈ λ(H sim ) or Z H ≈ cZ Hsim for some known constant c. The goal of such an encoding is to map the original Hamiltonian H onto a simulator Hamiltonian H sim where H sim is stoquastic and its low-energy dynamics are effectively captured by H. The upshot is that the encoded Hamiltonian H sim is directly computationally useful. If for a class of Hamiltonians H, simulator Hamiltonians H sim can be found, then the groundstate energy problem or the adiabatic computation of H is amenable to quantum Monte Carlo techniques and its complexity is in StoqMA. Additionally, the definition of encoding employed in [15] ensures that the partition function of the simulator Hamiltonian can be used to estimate that of the original Hamiltonian. An example is the mapping of any k-local Hamiltonian onto a real k + 1-local Hamiltonian [15].
Surprisingly, there are no known no-go's for the use of perturbative gadgets to map seemingly non-stoquastic Hamiltonians onto stoquastic simulator Hamiltonians. Even though mapping all Hamiltonians onto stoquastic Hamiltonians seems to be excluded from a computational complexity perspective, it is possible that there are Hamiltonians which are not stoquastic by a local basis change but for which a stoquastic simulator Hamiltonian exists.

Definitions
To make contact with some of the mathematics literature, we will use some standard terminology:

Definition 1 (Symmetric Z-matrix). A matrix is a symmetric Z-matrix if all of its matrix elements are real, it is symmetric, and its off-diagonal elements are non-positive.
Note that for a symmetric Z-matrix H the matrix exponential exp(−H/kT ) is entrywise non-negative for all kT ≥ 0 1 .
We reproduce the definition of a Hamiltonian being termwise stoquastic from [1] 2 . We will refer to it as stoquastic for simplicity. Given the discussion in the introduction and following [16], we also introduce the following definition (loosely stated): Another aspect of stoquasticity is the complexity of determining whether H is (computationally) stoquastic. The question is only clearly formulated when we restrict ourselves to specific mappings. For example, one can ask: is there an efficient algorithm for determining whether a 2-local Hamiltonian on n qubits is stoquastic in a basis obtained by single-qubit Clifford rotations. In this light, it was shown in [16] that for 3-local qubit Hamiltonians determining whether singlequbit Clifford rotations can indeed "cure the sign problem" is NP-complete. For two-local Hamiltonians the question how hard it is to decide whether a Hamiltonian is stoquastic by local basis changes is entirely open. The next section reviews our main results on this matter.

Summary of Results
In Section 2 we consider, as a warm-up, the two-qubit Hamiltonian problem and analyze under what conditions such a Hamiltonian is stoquastic by a local basis change. This turns out to be surprisingly complex. In Theorem 7 we show that one can determine whether a 2-qubit Hamiltonian is real under local basis changes by inspecting a subset of non-local invariants. In Section 2.1 we show that basis changes beyond single-qubit Clifford rotations are strictly stronger than Clifford rotations. For the general two-qubit case we show that a two-qubit Hamiltonian H is stoquastic iff one can find a solution to a set of degree-4 polynomials in two variables (the polynomials are quadratic in each variable), with details in Appendix C.
In Section 3 we move to the problem of determining stoquasticity by local basis changes for n-qubit 2local Hamiltonians. In general, the problem of deciding whether a multi-qubit Hamiltonian is stoquastic by a local basis change can be hard, since the local rotations bringing the two-qubit terms to symmetric Z-matrix form have to be mutually consistent. In addition, there is a subtlety involving the interplay of single and twoqubit terms. Definition 2 requires one to find a decomposition into 2-qubit terms D i each of which is a symmetric Z-matrix, i.e. 1-local terms can be allocated to different 2-qubit terms, leading to different decompositions D i . For a fixed basis, we show at the beginning of Section 3 that an efficient "parsimonious" strategy can find the optimal decomposition D i given a fixed basis (a similar idea is discussed in the Suppl. Material of [16]). If we allow basis changes, rotating the single-terms, the problem becomes potentially more complex. Proposition 8 shows that for an n-qubit 2-local Hamiltonian on a bipartite graph with identical terms everywhere, stoquasticity of the two-qubit case gives a sufficient condition for stoquasticity of the full Hamiltonian, and a necessary condition when restricted to local basis changes which are identical on each bipartition. Such classes of Hamiltonians are well motivated physically.
Our most important result is that for a subclass of nqubit Hamiltonians without 1-local terms, namely XYZ Heisenberg models on arbitrary graphs, we can find an efficient algorithm to determine stoquasticity by local basis changes. The Hamiltonian of the XYZ Heisenberg model for n qubits is Here the coefficients a uv P P , which can be different for each uv, are specified with some k bits. The theorem that we prove in Section 4 is: There is an efficient algorithm that runs in time O(n 3 ) to decide whether an n-qubit XYZ Heisenberg model Hamiltonian is stoquastic. The algorithm finds the local basis change such that the resulting H is a symmetric Z-matrix or decides that it does not exist.
This result relies partially on Lemma 9 (proven in Appendix D) which shows that when an XYZ Heisenberg model is stoquastic, it is so by a basis change induced by single-qubit Clifford rotations.
Going beyond the XYZ Heisenberg Hamiltonian, one can ask about the complexity of deciding stoquasticity by local Clifford rotations. Even though this is not the most general class of local basis changes, see for example Equation (7) in Section 2, our algorithm for the XYZ Hamiltonian suggests progress could most readily be made in this direction. We are able to prove that determining whether a Hamiltonian is real by local Clifford rotations is NP-complete: Note that this result does not imply that deciding stoquasticity by single qubit Cliffords or even general local basis changes for 2-local Hamiltonians is NP-complete. We give the proof of Theorem 5 in Appendix E since it is not central to the paper.

Warm-Up: Two Qubit Hamiltonians
Consider a general two qubit Hamiltonian where we may take a II = 0 without loss of generality. It is simple to show that Visually, the positive off-diagonal contributions of XX and YY terms need to be removed leading to a XX ≤ −|a Y Y |. Similarly, the positive off-diagonal contributions of IX and ZX terms need to be removed, leading to a IX ≤ −|a ZX | (and similarly a XI ≤ −|a XZ |).
The set of traceless symmetric Z-matrices forms a polyhedral cone C 2 , that is, a cone supported by a finite number of hyperplanes. In the case of 4 × 4 matrices, each matrix is determined by 9 parameters (three of them real and 6 of them nonnegative). Viewed as a polyhedral cone, C 2 has 12 extremal vectors [17] Appendix A discusses some of this structure and its generalization to two qudits.
Two-qubit Hamiltonians are conveniently parametrized by a 3 × 3 real matrix and two threedimensional real vectors: This implies that for 2-local qubit Hamiltonians it suffices to consider pairs of local SO ( To characterize the set of stoquastic 2-qubit Hamiltonians one needs to consider the orbit of C 2 under the local rotations. For a trace-zero Hermitian matrix acting on two qubits, there is a complete set of 18 polynomial invariants {I l } specifying the Hermitian matrix up to local unitary rotations, which we will refer to as the Makhlin invariants [18]. Note that the fact that there are 18 polynomial invariants can be consistent with the existence of only d inv = 15 − 6 = 9 independent nonlocal parameters in a two-qubit traceless Hermitian matrix [19,20]. One can see these 18 Makhlin invariants simply as a convenient set of invariants.
Stoquastic two-qubit Hamiltonians thus form a volume in the d inv = 9-dimensional non-local parameter space: the question is how the invariants are constrained in this stoquastic subspace.
We have not been able to express the stoquasticity of a two-qubit Hamiltonian in terms of inequalities involving the Makhlin or other invariants. However, we are able to capture the condition that H can be made real by local rotations in terms of the Makhlin invariants, see Theorem 7.
The convenience of the Makhlin invariants is that they can be expressed in terms of inner products and triple products constructed from β, S and P which are invariant under SO(3) rotations. Of particular interest are the triple product Makhlin invariants I 10 , I 11 and I 15 − I 18 . These take the following form (using the notation that ( a, b, c) stands for the scalar triple product a · ( b × c)): One can prove that Table 6 are equal to zero.

Theorem 7. A two-qubit Hamiltonian H is real under local unitary rotations if and only if all of the triple product invariants given in
The proof of this Theorem is simple in one direction, namely when H is real, then the triple product invariants of the corresponding β, S, P are zero. This can be seen by observing that realness implies that the vectors S and P lie in a common two-dimensional subspace and the repeated multiplication by β or β T keeps them in this subspace. Any triple of vectors, all lying in a twodimensional space have zero triple product.
In order to prove the other direction, one can observe To prove that zero triple product invariants implies realness under local rotations is more work since β can be degenerate and have rank less than 3, so the proof has to incorporate the various cases. We have included a full proof of Theorem 7 in the Appendix B.

Stoquastic Characterization: Beyond Clifford Rotations
The goal of this section is to better understand what freedoms to exploit when transforming a two-qubit Hamiltonian into a symmetric Z-matrix by local basis changes. The single-qubit Clifford rotations play a special role as local basis changes. The action of singlequbit Clifford rotations on the Pauli-matrices form a discrete subgroup of SO(3) representing the symmetries of the cube [21]. The Clifford rotations realize any permutation of the Paulis (the group S 3 ) as well as signflips σ i → ±σ i (with determinant 1). However, when deciding if a 2-qubit Hamiltonian is stoquastic, it does not suffice to consider only the local Cliffords.
We can consider a concrete example It is easy to argue that no single-qubit Clifford basis change can make this Hamiltonian into a symmetric Zmatrix. Since H has to be real, no permutations to Y -components are allowed. A sign-flip which makes the XX term negative also messes up the sign of one of the single-qubit X terms. Interchanging X and Z for both qubits leads to the same problem for the ZZ and Z terms. Applying a sign flip on the Z of one of the qubits, followed by interchanging X and Z on that same qubit leads to the requirement that −a Z ≤ −1 and −a X ≤ −|a XX |. This is not satisfiable as long as 0 < a Z < 1.
However, we could create off-diagonal XZ and ZX terms by a single-qubit non-Clifford rotations. For example, on both qubits one can apply a π-rotation around the Z-axis followed by a π/4-rotation around the Y -axis mapping For small a XX 1, and sufficiently large |a X − a Z |, these inequalities can certainly be satisfied.
In Figure 1 we show in what region of the parameter space (a X , a Z , a XX ) H is stoquastic using the general set of inequalities given in Appendix C.
Since local Cliffords are insufficient, we need a more general strategy. A general analytical strategy to determine if a two qubit Hamitonian is stoquastic, and in what basis, is as follows. We assume that one first applies local rotations to the two-qubit Hamiltonian H to bring it to a standard form where β is diagonal: A trivial consequence of the previous section is that H is real under local unitary rotations if and only if the diagonalizing transformation lead to the 1-local terms being real: Note that when β has a degeneracy, for example when |a XX | = |a Y Y |, an additional rotation may need to be performed on the degenerate space to retrieve the above form.
Without loss of generality we can assume that a ZZ ≥ a XX ≥ 0. This is because we can always perform sign flips on the eigenvalues of β by Pauli rotations, and dump any extra sign into a Y Y so as to preserve the determinant.
As a simple special case, when S = P = 0, H can always be made into a symmetric Z-matrix by Clifford rotations, i.e. choosing an appropriate choice of permutations and sign-flips on β such that a XX ≤ −|a Y Y |. So let us assume that at least one of the vectors S and P is not zero. Furthermore, if β = 0 then H is stoquastic, since S and P may be freely rotated. As a second special case, if a ZZ = a XX = 0 and a Y Y = 0 then any transformation of our Hamiltonian into a symmetric Zmatrix will have to make a Y Y = 0. Therefore either the Hamiltonian is not stoquastic or a permutation can be applied preserving realness and setting a Y Y = 0 and a ZZ = 0. Thus we can assume a ZZ > 0 and normalize our Hamiltonian such that a ZZ = 1.
We wish to know then if there exists a pair of SO(3) where β , S , P are associated with a symmetric Z-matrix Hamiltonian. Any such transformation must preserve the realness of the Hamiltonian. If at least one of the vectors S and P is rank-2 (if this is not the case, see the extra step in Appendix C) then the only realness-preserving transformations which may be performed are SO(3) rotations in the X-Z subspace combined with reflections in the X-or Z-axis, given by This gives a β , S , P . Recalling Proposition 6, three inequalities must be satisfied in order for the rotated H with β , S , P to be a symmetric Z-matrix. As we show in Appendix C these inequalities can be re-expressed as systems of two-variable polynomials which are at most quadratic in either variable, so analytic solutions to their roots can be constructed and solutions can be found using graphical methods. However, the complexity of these inequalities makes it unclear whether their interest goes beyond that of numerically finding a set of local basis changes.

The Stoquastic Problem for 2-local Hamiltonians
Consider a general n-qubit 2-local Hamiltonian where u and v are the vertices of the interaction graph corresponds to a qubit so that H v are the 1-local terms (corresponding to S and P in the two-qubit case). Each edge e = uv ∈ E is represented as H uv and can be associated with the 3 × 3 matrix β uv for that edge. We can define the set of 2-local symmetric Z-matrices as the polyhedral cone C n which is generated by the extremal vectors (viewed in a n-qubit Hilbert-Schmidt space, with Is appended) of each of the 2-qubit cones C 2 . The cone C n is thus the set of n-qubit 2-local Hamiltonians which are stoquastic with respect to Definition 2 without further basis changes. If we do not allow for basis changes, it is simple to determine whether a n-qubit Hamiltonian is real by checking the realness conditon.
So for the next step we assume that H is real and want to decide whether H ∈ C n .
It is important to note that for any v, C uv 2 and C vu 2 intersect: for any u, u = v, the extremal vectors obey This non-empty intersection shows the freedom in the decomposition D i for a fixed basis. The parsimonious strategy below shows how to find a decomposition D i of a n-qubit Hamiltonian H in terms of 2-local terms each of which is a 4 × 4 symmetric Z-matrix or decide that it does not exist.
Let H X u be the 1-local term proportional to X u and similarly H Z u . First, note that a real H can only be in C n when for all u, H X u ≤ 0. This follows directly from Proposition 6 demanding that a IX , a XI ≤ 0. Hence we assume this to be the case (otherwise decide not stoquastic in the given basis).
Efficient Parsimonious Strategy: Repeat the following for all edges uv in H.

Given the current Hamiltonian H, pick a pair of vertices u and v and consider h(α, β)
which includes all current singlequbit terms which act on vertices u and v and wlog In the latter case, decide that H ∈ C n and exit. Note that when h(α, β) / ∈ C 2 then ∀α ≤ α, β ≤ β, h(α , β ) / ∈ C 2 since a IX ≤ −|a ZX | and a XI ≤ −|a XZ |, aka it is easier to satisfy the stoquastic condition for large α, β, more negative 1-local Xterms. Hence the goal is to use as little α and β as possible, be parsimonious, so that a lot of −X u and −X v will be left over for another edge. This is then clearly an optimal strategy. Let us next consider the problem of admitting local basis changes. A particular simple example for which a resolution of the two-qubit problem analyzed in the previous section suffices, is the following: Furthermore, H = uv∈E h uv where h uv acts with both one and two local terms on sites u and v, and

If the two-qubit Hamiltonian h is stoquastic, then H is stoquastic. If h is not stoquastic then H will not be stoquastic under any local basis change which acts identically on all qubits in a partition.
Bipartite graphs include linear arrays, square lattices, cubic lattices and hexagonal lattices, all of which are very natural structures to consider.
Suppose h were not stoquastic, but there existed a pair of unitaries U A and U B such that a decompo- cannot be a symmetric Zmatrix, and so D uv = h uv . However both D uv and h uv must share the same purely two local parts. As such, D uv must differ from h uv by its 1-local terms, and in particular one or both of the −X u and −X v terms. Since h uv is not stoquastic, D uv must have more support on either −X u or −X v than h uv . Suppose wlog that D uv has more support on −X u , then, by the parsimonious reasoning outlined earlier, there must exist some D ux which has less support on −X u than h ux does, and consequently D ux is not a symmetric Zmatrix, leading to a contradiction.

Efficient Algorithm for the XYZ Heisenberg Models
In this section we will outline and detail the proof of our main Theorem, Theorem 4. Consider an XYZ Heisenberg model Hamiltonian on n qubits and arbitrary interaction graph G: with For each H uv , S = P = 0 and The aim is to show an efficient algorithm (which we will show has O(n 3 ) steps) which either decides that H is not stoquastic by local basis changes or finds it is stoquastic (and also gives you the local basis change).
As a result of S and P being zero for every term H uv , it follows from Proposition 6 that H uv will not be a symmetric Z-matrix unless β uv is in diagonal form. Thus we can say that H is stoquastic if and only if there exists a set of SO (3) which is equivalent to We emphasize that we order the indices Z → 1, X → 2, Y → 3 to denote the entries ofβ. The proof of Theorem 4 now has the following outline. First, we prove a Lemma which establishes that it suffices to consider single-qubit unitary Clifford rotations in order to determine whether H is stoquastic. The action of these single-qubit Clifford unitaries on a matrix β uv is simply to permute the diagonal elements and add ± signs to the diagonal elements, which is the only freedom when β is to remain diagonal. Since we know that in order for H to be stoquastic, β should remain diagonal, it seems natural that the set of Clifford rotations indeed suffices, but the proof which can be found in Appendix D requires some actual work.

Lemma 9 (Single-Qubit Clifford Rotations Suffice). The XYZ Heisenberg Hamiltonian is stoquastic by local basis changes if and only if it is stoquastic by single-qubit Clifford rotations. To be precise: given a set of diagonal
for all u and v. Before we start with proving the main theorem, we observe the fact that a solution set with Π u ∈ O(3) can be always be converted into a solution set for which all Π u ∈ SO(3). More precisely, for any solution set ) and Π u ∈ S 3 , which obeys conditions (12) and (13), where some 11 is irrelevant. By Lemma 9, in order to determine whether H is stoquastic, one needs to consider whether there is a set of permutations {Π u } such thatβ uv = Π u β uv Π T v obeys Equation (14). If such set of permutations cannot be found, one decides that H is not stoquastic. If such sets of permutations can be found, one needs to determine whether for any of these solution sets of permutations one can find sign flips {R u } such that Equation (15) is also obeyed.
In order to determine the permutation solution set {Π u }, we first establish in a simple lemma (Lemma 11) that for any edge β uv which has rank 2 or 3, -we can call these rigid edges-, any solution set needs to have Π u = Π v . Hence for each rigid connected component (see Definition 10) one has only 6 possible permutations Π ∈ S 3 to choose from and we can iterate over these choices to determine a candidate permutation solution set, see Proposition 12. For the remaining rank-1 edges, which can be labelled with colors 1, 2 or 3 based on which diagonal element of β uv is non-zero, Condition (14) implies that none of these edges can haveβ uv 33 = 0, hence the permutations have to eliminate edges of color 3 while keeping β diagonal and being members of the candidate permutation subset. By eliminating the freedom residing in clusters of internal vertices which only connect to vertices by single-color edges (see Definition 13), the problem of determining the permutation solution set is then shown to map onto the problem of finding the satisfying solutions of some (disjoint) Ising problems (Lemma 14). In Section 4.1 we then go back to the determining {R u } on the basis of these permutation solutions and argue that the additional freedom in internal vertices does not cause any complexity blow-up.
First, we define rigid connected components and the flexible quotient multi-graph: Definition 10 (Rigid Connected Components and Flexible Quotient Multi-Graph of G). Take the graph G = (V, E) associated with a Hamiltonian H and remove all edges associated with β uv of rank 1. The resulting graph G − is empty or has connected components Γ − i , i = 1, . . . , K. To each connected component Γ − i we now add the rank-1 edges in G connecting vertices in Γ − i : these graphs will be called Γ i , i = 1, . . . , K and they are the rigid connected components of G. We can define the flexible quotient multi-graph G Q = (V Q , E Q ), see Fig. 2, with a vertex i ∈ V Q for each rigid connected component Γ i (thus |V Q | = K). For each rank-1 β uv edge uv ∈ E with u ∈ Γ i and v ∈ Γ j with j = i, the graph G Q has an edge between i ∈ V Q and j ∈ V Q . Note that G Q can have multi-edges.
Lemma 11 (Permutations are identical on rigid connected components). Given a rigid connected component Γ and a solution to that graph {Π u = R u Π u }, there exists a permutation Π such that Π u = Π for every vertex u in Γ. Figure 2: An illustration of the mapping of a graph G to its flexible quotient multi-graph GQ, with red edges corresponding to rigid edges and black edges associated with rank-1 β uv .
. Since β uv mm = 0 for two of the three values of m, Π u = Π v . This holds for every pair of connected vertices in the connected components, hence for all u ∈ Γ, Π u = Π for some Π.
We now work with the flexible quotient multi-graph G Q . If this graph G Q has multiple connected components, we treat each component separately, so assume it has a single connected component.

Definition 12 (Candidate Permutation Set).
For each rigid connected component Γ i (associated with a vertex of G Q ), we consider all 6 permutations Π ∈ S 3 on its vertices and determine the subset S i ∈ S 3 for which all edges β uv satisfy Condition (14). If any S i is the empty set, we decide that H is not stoquastic. This iteration over at most 6n 2 permutations is clearly efficient. We call {S i } i=K the candidate permutation solution set. If there are no rigid connected components in G, the candidate solution set is the set of all permutations S ×n 3 .
Let G Q thus be a multi-graph with an associated candidate permutation set {S i } K i=1 such that each simple edge or loop has one of three different indices or colors x = 1, 2, 3. This coloring is determined by which diagonal element in β is non-zero. The problem of determining whether there exists a set of permutations Π i , one at each vertex i, such that for each simple edge (i, j) of the multi-graph, one has is equivalent to determining whether there exists a set of permutations Π u which obey Condition (14) in the original graph G. Conditions (16) and (17) ensure that every edge will be permuted into an edge not labelled by 3 which is required due to each edge being of rank 1 and Condition (14). Condition (18) ensures that every edge is mapped to an edge which is still in diagonal form (Condition (12)) and the last Condition (19) verifies that the permutations are consistent with the internal constraints in a connected rigid subgraph. The next step is to alter the multi-graph G Q to an effective multigraph G eff Q that will be analyzed in Lemma 14. Definition 13 (From G Q to G eff Q : Removing Single-Color Clusters). Given the flexible multi-graph G Q associated with a graph G and a Hamiltonian H. If a vertex in G Q is incident to edges of three different colors, 1,2, 3, then H is not stoquastic due to constraints (16) and (17). If not, construct the effective multi-graph G eff Q as follows. For each color x = 1, 2, 3 define the subset of internal vertices V x ⊆ V as those which are only incident to edges of color x. V x can be composed of dis- Note that each permutation solution set of G eff Q uniquely determines the map x → x for each cluster V c x due to Condition (18). This also means that the permutations on these internal vertices are not necessarily fixed, for each internal vertex j there are two permutations Π 1 j (x) = x and Π 2 j (x) = x which work as possible permutations as long as Π 1 j , Π 2 j ∈ S j . Since these permutations have the same permutational action, it is irrelevant which one to take, but when we expand internal vertices in G Q to rigid connected components again, we will discuss the effect of these freedoms in efficiently solving for {R u } in section 4.1. Proof. Since every vertex is incident to edges with only two colors, there are three types of vertices, namely 12, 23 and 13 vertices, see Figure 3. Then, for each type of vertex, only two possible permutations obeying conditions (16)-(17) exist and we can label these choices as Ising spins ±1, as illustrated by each pair of antipodal points in Figure 3. The last condition, Condition (19), which expresses the candidate solution set, imposes a restriction on these Ising spins, i.e. they can fix the Ising spin to be +1 or −1, modeled as some h i S i = −1. Furthermore, note that if we choose one of the two permutations for one vertex i which is connected by an edge x to a vertex j, then this uniquely fixes the permutation at vertex j due to Equation (18). So for each edge, one can define a ferromagnetic (J ij = 1) and anti-ferromagnetic Ising interaction (J ij = −1) as follows:

Lemma 14 (Ising Problem). Let G eff Q be a multi-graph on n vertices such that each vertex is incident to edges of precisely two colors and a candidate permutation solution set for each vertex. The problem of determining whether there are permutations of the colors on the vertices of G eff Q which obey conditions (16)-(19) is equivalent to the question whether there is a {S
• Any color-3 edge with one incident 13-vertex and one incident 23-vertex will act as an antiferromagnetic bond J ij = −1 on the {+1, −1} labels of the permutations in order satisfy conditions (16)- (18). This can be seen easily in Figure 3 by noting the labelling of the vertices connected to color-3 edges.
• Any color-3 edge where both incident vertices are either 13 or 23 will act as a ferromagnetic bond J ij = +1.
• Any color-2 or color-1 edge, regardless of its incident vertices will act as a ferromagnetic bond J ij = +1.
Once a single permutation is fixed, the Ising interactions determine all other permutations, and one can verify whether these permutations are consistent with the restriction due to the candidate solution set. Since a single vertex only admits at most two possible assignments of permutations, this means that one only needs to make at most two of these checks.

Determining The Sign Flips {R u }
If there is a stoquastic basis for H, then Lemma 14 will give either one or two set of permutations {Π j } for G eff Q . These solution sets can then be extended to possible permutations on G Q which include the internal vertices in clusters. It may be the case, depending on the candidate permutation sets, that for some of the internal vertices j in a cluster V c x , two possible permutations Π j (x) = x and Π j (x) = x are valid permutations  (16), (17) and (18) are satisfied. The schematic shows that a permutation on a vertex which is incident to two different edges x = y can only be one of two choices (antipodal points in the figure) which can be labelled as ±1. In addition, once we fix a permutation for a vertex, only one choice of permutation is allowed for its adjacent vertex: for example if we take a 12-vertex and fix its permutation to be (1)(2)(3), any neighboring 12-vertex must also be assigned (1)(2)(3) while any neighboring 23-vertex must be assigned (13)(2).
to make H stoquastic. We call these permutations solutions on internal vertices internal permutations. In addition, the action of one of these permutations on a vertex will extend to the rigid connected component associated with that vertex, leading to some possible permutation solution sets {Π u } |V | u=1 for the whole graph G.
Let's first establish that given a fixed permutation solution set {Π u }, it is simple to verify whether one can satisfy Condition (13), i.e. do there exist sign-flips {R u } such that (R u Π u β uv Π T u R T u ) 22 ≤ 0. For a given set {Π u }, we construct a new graph G 2 from G. G 2 has the same vertices as G but we keep only the edges for which (Π u β uv Π T u ) 22 = 0. With each edge we associate the sign of [Π u β uv Π T v ] 22 as a ±1 sign. Note that G 2 is not necessarily connected. 22 > 0 then in order to satisfy Condition (15) we want to apply a sign flip R = diag(1, −1, 1) to either vertex u or vertex v. If [Π u β uv Π T v ] 22 < 0 we want to either apply the sign flip to both vertices u and v, or not apply a sign flip to either. Again, this is precisely a question of satisfying all Ising interactions in a classical Ising problem. Here negative edges are mapped onto an anti-ferromagnetic bond, and positive edges are mapped onto a ferromagnetic bond. If G 2 has several disconnected components, then each component will have two possible independent sign solutions which can be efficiently enumerated.
If our graph G Q has clusters with internal vertices, on which there can be more than one internal permutation solution, one has to be careful to argue that one does not end up with an exponential number of such Ising problems due to the multiplicative number of choices of permutation assignment for the whole set of internal vertices. In order to argue that this is not the case, we need the following observations: • If a rigid connected component Γ i ∈ G is represented by an internal vertex i in a cluster V c x where x is mapped to 1 by the permutation solution of G eff Q , then in the graph G 2 all the vertices of Γ i are disconnected from the other vertices in the G 2 graph. This implies that each Γ i in G representing an internal vertex i in V c x corresponds to a disconnected component in G 2 . For such a disconnected component one has to consider whether one can find a set of sign-flips for at most two internal permutations, so this requires checking at most two Ising problems, and then disregarding that disconnected component while checking the rest of the internal vertices.
• A rigid connected component cannot be a vertex in a cluster V c x where x is mapped to 3 by the permutation solution of G eff Q since the permutation solution has to eliminate simple edges of color 3 by conditions (16) and (17). Hence this case is excluded. Thus, despite the large number of possible permutation assignments to the full set of internal vertices, the impact of a given choice of permutation for any individual internal vertex i on the subsequent problem of assigning sign-flips is either non-existent, or isolated to only those vertices in G associated with the rigid subgraph Γ i , and can be considered independently of the other choices. Therefore the number of Ising problems one will need to consider is bounded. To summarize, the algorithm requires the following steps, with the following cost, where n is the number of vertices in the graph G.
(1) First the rigid connected components are identified. A candidate permutation set is decided for each rigid connected component The flexible quotient graph is constructed by modding out the rigid connected components. Each edge is coloured according to what entries are non-zero in their β matrix.   (5) Solving the Ising model gives two possible assignments of signs. Some of these assignments may be incompatible with candidate permutation sets of a given rigid connected component. (6) With a sign assignment, we may use figure 3 to assign permutations.
This fixes all permutations in G Q except for isolated vertices.   • Identifying rigid connected components: O(n).
• Constructing candidate permutation sets for the rigid connected components: O(n 2 ).
• Constructing G eff Q by first identifying and removing internal vertices, which takes O(n 2 ) time, then linking up boundary vertices, which also takes O(n 2 ) time.
• Solving the Ising problem for G eff Q : O(n 2 ). • The above Ising problem produces at most two possible sets of permutation assignments on G eff Q , each of which correspond to multiple possible sets of distinct permutation assignments on G. However this number is at most 2 times the number of internal vertices in V x for which x is mapped to 1. So the number of possible permutation assignments on G is O(n). Hence we conclude that the problem of determining whether there is a basis change {D u } such that H is stoquastic can be done efficiently in worst case O(n 3 ) time. It should be noted that a more careful analysis would probably find a tighter worst-case of O(n 2 ), since the worst-case of O(n 2 ) for deciding {R u } assumes all vertices in G 2 are connected, while the worst-case of O(n) for the number of permutation assignments would imply that almost none of the vertices in G 2 are connected. For the sake of clarity, we have included an illustrated overview of the algorithm in tables 1 and 2.

Discussion
There remain many open questions with regards to stoquasticity. It would be interesting to capture the twoqubit (or two-qudit) stoquastic space in terms of inequalities involving invariants. It would be desirable if the problem of deciding whether a Hamiltonian is computationally stoquastic is easier when allowing more transformations than just local basis changes. Such simplification would also be welcome for the two-qubit case.
Various of the proofs in this paper are tedious and involved since the reasoning is different for different subcase Hamiltonians. We did not see a way of simplifying these arguments. It is an interesting question whether one could supply such proofs using a proof assistant, automatically ensuring their validity.
The presence of single-qubit terms causes our algorithm for deciding the stoquasticity of the XYZ Heisenberg model to fail. Furthermore, the existence of singlequbit terms makes the class of single-qubit Clifford rotations strictly weaker than that of single-qubit rotations. Nevertheless, it seems likely that a generalization of the algorithm described here could be constructed to decide stoquasticity under single-qubit Clifford rotations of an XYZ model with single qubit terms included. It remains an open question whether it is computationally hard to decide whether such a Hamiltonian is stoquastic by local basis changes beyond single-qubit Clifford rotations.
We expect that the problem of deciding whether a Hamiltonian is stoquastic by local basis changes is easy for n-qubit Hamiltonians with an underlying line or tree interaction graph G by employing a combination of an efficient dynamic programming strategy with the parsimonious strategy (i.e. discretize the set of local basis changes, show how the optimization of a partial problem giving a basis change at the boundary can be used to solve a next-larger partial problem etc.).
An intriguing idea is that results showing hardness of deciding whether a Hamiltonian is stoquastic could be used for quantum computer verification. A general adiabatic computation with frustration-free Hamiltonians requires a quantum computer to run it, but let's assume that Alice, the verifier, picks a stoquastic frustrationfree Hamiltonian, but then hides this stoquasticity by local basis changes or, say, the application of a constant depth circuit. Bob with a quantum computer can run the adiabatic computation and provide Alice with samples from its output which she can efficiently verify to be statistically correct using the algorithm in [8]. On the other hand, a QC-pretender Charlie may have to either find the constant-depth hiding circuit or classically simulate the quantum adiabatic algorithm, and either task could be too daunting.

Funding
We acknowledge support through the EU via the ERC GRANT EQEC No. 682726. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation.

Acknowledgements
BMT would like to thank David DiVincenzo for interesting discussions. BMT and JDK would also like to thank Sergey Bravyi as well as Christophe Vuillot for helpful discussions.
It is understood that the set of symmetric Z-matrices forms a polyhedral cone. For convenience and concreteness we present the structure of this cone for operators acting on C d 2 in terms of generalizations of the Pauli matrices.
Any Hermitian operator acting on C d ⊗ C d with basis elements {|i | i ∈ 0, ..d − 1}, can be written as : with a xy ij;mn ∈ R. Here the λ D,A,S ij matrices with labels S (symmetric), A (anti-symmetric) and D (diagonal) are a modified set of generalized Gell-Mann matrices [22] given by: with e i = |i i| − 1 d I. Note that Tr[λ x ij λ y mn ] = δ im δ jn δ xy d, so that the λ matrices can be thought of as spanning vectors of the Hilbert-Schmidt vector space. We can thus employ the vector inner product of two bipartite matrices A and B acting on C d ⊗ C d as: The necessary and sufficient conditions for H to be a symmetric Z-matrix (see Definition 1) are (1) ∀i, j, m, n im|H|jn ∈ R and (2) im|H|jn ≤ 0 when either i = j or m = n (or both).
The λ-basis is convenient because of the following properties. The diagonal part of H only has support on vectors of the form λ D ⊗ λ D while the off-diagonal part has no contribution from λ D ⊗ λ D . Furthermore, the imaginary part of H only has support on vectors of the form λ A ⊗λ S , λ A ⊗λ D , λ S ⊗λ A and λ D ⊗λ A , while both the real and the diagonal part of H have no contribution from these vectors. Thus the vector H has support on three orthogonal subspaces, the diagonal subspace D spanned by vectors λ D ⊗ λ D , the imaginary subspace I, spanned by vectors λ A ⊗ λ S , λ A ⊗ λ D , λ S ⊗ λ A and λ D ⊗ λ A , and the off-diagonal real subspace R spanned by the remaining vectors.
We may write H = H D +H R +H I where H D , H R , H I are the vectors in these three subspaces.
The first ("realness") condition thus corresponds to demanding that H I = 0. By the linear independence of the λ-basis it implies that ∀i, j, m, n, a AS ij;mn = a SA ij;mn = a AD ij;mn = a DA ij;mn = 0. (23) The second ("negativity") condition then says that for all i = j or m = n (or both): This condition thus specifies some of the facets of the symmetric Z-matrix cone. When a vector |im jn| is interpreted as the normal vector of an oriented hyperplane in Hilbert-Schmidt space, the negativity condition can be thought of as the requirement that H R lies within the union of half spaces defined by these oriented hyperplanes.
To write down the inequalities describing these facets, we can use that One can then write down three types of inequalities: • case 1: i = j and m = n There is an interesting geometric fact here. Consider for instance Inequality (25). This is illustrated in figure 4 for d = 2, 3. In the case of two qubits λ S 01 = X, λ A 01 = Y , λ D 01 = Z, e 0 = Z/2 and e 1 = −Z/2. For the realness condition we require that a IY = a Y I = a Y Z = a ZY = a Y X = a XY = 0. For the negativity condition we retrieve the inequalities: Together these (in)equalities are given in condensed form in Proposition 6.

B Proof of Theorem 7
Recall the definitions of When the triple product invariants are zero, Γ L and Γ R each contain coplanar vectors. By Proposition 15 (see below) we know that S is in the span of at most two left singular vectors of β and P is in the span of at most two right singular vectors of β.
If S, ββ T S and [ββ T ] 2 S are linearly independent, then βP must be in the same span as S. If they are not linearly independent then S is a left singular vector of β, with unit vectorê S . It follows that βP ∈ span(ê S ,ê S ⊥ ), and ββ T βP ∈ span(ê S ,ê S ⊥ ). This implies that ββ Tê S ⊥ ∈ span(ê S ,ê S ⊥ ). Sinceê T S ββ Tê S ⊥ = 0 it follows that ββ Tê S ⊥ ∈ span(ê S ⊥ ) andê S ⊥ is a left singular vector of β. Therefore there always exists a pair of left singular vectors of β such that S and βP are in their span. By Proposition 16 (see below) this implies that S is spanned by at most two left singular vectors of β and P is spanned by the corresponding right singular vectors. By Proposition 17 (see below) we see that H is real under local unitary rotations. Proof. Let S ∈ span(ê L x ,ê L y ). Assume βP is in the same span as S, then β T βP, [β T β] 2 P ∈ span(ê R x ,ê R y ). If β T βP and [β T β] 2 P are not proportional to one another then P ∈ span(β TêL x , β TêL y ). If β T βP and [β T β] 2 P are proportional to one another then they are proportional to a right singular vector of β,ê R P . If β is full rank then β T β is invertible and P ∝ê R P ∈ span(ê R x ,ê R y ). If β is not full rank and P ∝ê R P then P ∈ span(ê R P ,ê R 0 ) with βê R 0 = 0. Consider the left singular vector expansion of S into an orthonormal basis: We require the following three vectors to be coplanar  It is easy to see that S and P are spanned by at most two left singular vectors of β That concludes the proof for one direction of the biconditional.
Suppose S and P are spanned by two or fewer left and right (respectively) singular vectors of β which share the same singular values. Then S and P are expressible as

C Polynomial Inequalities for Stoquasticity of a 2-qubit Hamiltonian
The rotated two-qubit Hamiltonian H is described by Three inequalities must be satisfied in order for H to be a symmetric Z-matrix: Note that the parameters γ 1 and γ 2 have fallen out of these inequalities.
If the above inequalities can be satisfied, then the Hamiltonian is stoquastic. If they cannot be satisfied then the Hamiltonian is not stoquastic, barring a single case as follows. If either {a XI , a IX } = {0, 0} or {a ZI , a IZ } = {0, 0}, then an additional test of the same form must be applied after first performing a permutation such that the pair of zero coefficients are now in the Y position. For example, if {a XI , a IX } = {0, 0} then one must also perform the above test on the diagonal form: The inequalities (30), (31) and (32) can be broken up into six inequalities in order to ignore the absolute values: We wish to map these trigonometric inequalities into polynomial inequalities of two variables by defining x 1 = tan(θ L ) and x 2 = tan(θ R ). However the mapping will depend on the value of cos(θ 1 ) and cos(θ 2 ), and so breaks up into several cases. Case 1: cos(θ L ), cos(θ R ) = 0 Let δ L = Sign(cos(θ L )) and δ R = Sign(cos(θ R )) then Case 2: cos(θ L ) = 0, cos(θ R ) = 0 Let δ L = Sign(sin(θ L )) and δ R = Sign(cos(θ R )) then 2 )a 2 ZI ≥ 1 Case 3: cos(θ L ) = 0, cos(θ R ) = 0 Let δ L = Sign(cos(θ L )) and δ R = Sign(sin(θ R )) then Case 4: cos(θ L ) = 0, cos(θ R ) = 0 Let δ L = Sign(sin(θ L )) and δ R = Sign(sin(θ R )) then All of these inequalities are at most quadratic in any single variable x 1 or x 2 , so analytic solutions to their roots can be constructed and solutions can be found using graphical methods. If there exist values of x 1 and x 2 such that for some choice of δ L , δ R ∈ {−1, 1} at least one of the above cases is satisfied, then H is stoquastic.
If not, then we may say that H is not stoquastic save the exception described above.
This concludes the complete characterization of stoquasticity of a two-qubit Hamiltonian H.

D Proof of Lemma 9
Consider a set of 3 × 3 matrices {β uv } indexed by the pair uv which are diagonal in some fixed basis {ê 1 ,ê 2 ,ê 3 }. Here we establish that whenever there exists a set of SO (3) where by signed permutations we mean matrices of the form Π u = Π u R u , with Π u a permutation and R u a diagonal matrix with diagonal elements ±1. Therefore whenever considering the class of sets of SO(3) transformations which preserve the diagonality of the matrices {β uv }, it suffices to consider only the signed permutations. Of course this is trivially true when considering a single diagonal matrix, but is less obvious when considering a set of matrices.
The proof is achieved by introducing a consistent mapping from a given SO(3) matrix O to a signed permutation Π(O), which we call the Π-reduction. This mapping will have the property that for a diagonal We can associate each of the monomials in the right hand side of the above equation with a matrix which is zero for all terms in the matrix which do not contribute to that monomial, and sign(x) for all terms x which do contribute to that monomial. So for example: We may also impose a lexicographical ordering on the matrices, so that the matrix associated with bdi is lower in the ordering than the matrix associated with bf g. This is to ensure that the D-reduction has a consistent definition, and the particular choice of ordering is not important. We may then define the Π-reduction as follows: Note that the determinant of O is +1, so there always exists a non-zero monomial. Finally, note that the Π-reduction is always a signed permutation, but not always an SO (3)   Clearly, the realness-under-Clifford rotations problem is in NP, since a prover can give the verifier the singlequbit Clifford rotations which make the terms in H real.
To prove the problem to be NP-hard, we show that if one can efficiently solve the realness-under-Cliffords problem, then one can also efficiently solve an NPcomplete subclass of the exact cover problem, called restricted exact cover by three-sets.