# Automated detection of symmetry-protected subspaces in quantum simulations

Caleb Rotello<EMAIL_ADDRESS>Department of Physics and Quantum
Engineering Program, Colorado School of Mines, Golden CO, 80401  National
Renewable Energy Laboratory, Golden CO, 80401  Eric B. Jones
<EMAIL_ADDRESS><EMAIL_ADDRESS>Infleqtion, Louisville CO,
80027  National Renewable Energy Laboratory, Golden CO, 80401  Peter Graf
<EMAIL_ADDRESS>National Renewable Energy Laboratory, Golden CO, 80401
Eliot Kapit<EMAIL_ADDRESS>Department of Physics and Quantum Engineering
Program, Colorado School of Mines, Golden CO, 80401

(August 27, 2024)

###### Abstract

The analysis of symmetry in quantum systems is of utmost theoretical
importance, useful in a variety of applications and experimental settings, and
is difficult to accomplish in general. Symmetries imply conservation laws,
which partition Hilbert space into invariant subspaces of the time-evolution
operator, each of which is demarcated according to its conserved quantity. We
show that, starting from a chosen basis, any invariant, symmetry-protected
subspaces which are diagonal in that basis are discoverable using transitive
closure on graphs representing state-to-state transitions under $k$-local
unitary operations. Importantly, the discovery of these subspaces relies
neither upon the explicit identification of a symmetry operator or its
eigenvalues nor upon the construction of matrices of the full Hilbert space
dimension. We introduce two classical algorithms, which efficiently compute
and elucidate features of these subspaces. The first algorithm explores the
entire symmetry-protected subspace of an initial state in time complexity
linear to the size of the subspace by closing local basis state-to-basis state
transitions. The second algorithm determines, with bounded error, if a given
measurement outcome of a dynamically-generated state is within the symmetry-
protected subspace of the state in which the dynamical system is initialized.
We demonstrate the applicability of these algorithms by performing post-
selection on data generated from emulated noisy quantum simulations of three
different dynamical systems: the Heisenberg-XXX model and the $T_{6}$ and
$F_{4}$ quantum cellular automata. Due to their efficient computability and
indifference to identifying the underlying symmetry, these algorithms lend
themselves to the post-selection of quantum computer data, optimized classical
simulation of quantum systems, and the discovery of previously hidden
symmetries in quantum mechanical systems.

††preprint:

##  I Introduction

The analysis of symmetry is a central tool in physics and has enabled some of
the most profound discoveries in the field. Noether’s Theorem famously
connects the symmetries of a system’s action with conservation laws to which
that system’s equations of motion are subject [1]. Generally, the analysis of
symmetry, or the breaking thereof, allows one to constrain theories [2], solve
equations of motion more efficiently [3], and identify phases of matter [4].
Applications of symmetry analysis in quantum information include, but are not
limited to: quantum error correction [5], error mitigation on quantum hardware
[6, 7, 8, 9], and quantum machine learning model design [10].

Quantum computing can efficiently simulate quantum dynamics in regimes where
classical simulation becomes impossible [11]. However, current quantum
processors operate in a regime severely constrained by noise, with error rates
not yet sufficiently below most error correction thresholds [12]. Error
mitigation will therefore be critical in the interim before fault-tolerant
architectures can be scaled [13, 14]. In recent work, and despite its
limitations [15], the technique of post-selection has proven useful to
mitigate errors and extract useful results from quantum simulation experiments
(see e.g., [16]). Post-selection works by identifying measured states that
could only have come from error processes and excluding them from the
statistics used to calculate output quantities. The most obvious example is a
conserved quantity such as particle number. In such a case, any measured state
that does not preserve the conserved quantity must be the result of errors.
Due to the connection established by Noether, a more fundamental way to
describe post-selection is with respect to symmetry. According to this
description, post-selection works by checking the eigenvalue of the
simulation’s “fiducial” (i.e., initial) state under the symmetry operator
against the corresponding symmetry operator eigenvalues (e.g., value of a
conserved quantity) of individual measurement results in the dynamically-
generated output state of the simulation. If a particular measurement outcome
registers a different eigenvalue under the symmetry operator than the fiducial
state does, then the measurement is a result of error and can be discarded.
This procedure is restricted in the scope of its application, as symmetries of
a quantum system and their corresponding operators (i.e., conserved
quantities) are typically either engineered into the dynamics “by hand” or
identified by clever theoretical intuition. For a generic quantum system, the
relevant symmetry operator(s) may not be obvious a priori and may be difficult
to identify. Being able to perform post-selection in a manner that does not
require explicit identification of a symmetry operator would greatly increase
the technique’s applicability; this is the subject of this paper.

As a corollary, such an operator-free method for error detection also enables
additional applications such as more efficient classical simulation of quantum
systems via computational basis state reduction. For example, particle number
conservation in hardcore boson models can be used both for post-selection in
quantum simulation and for reducing the basis state set size from $2^{n}$ to
$\binom{n}{N}$, where $n$ is the number of lattice sites and $N$ the number of
particles, in classical simulations [17]. Interestingly, the identification of
symmetry or conserved quantities in some instances can make classical
simulation so efficient that it can obviate the need for quantum computation
altogether [18]. Finally, in certain special cases, one may be able to infer
the explicit form of a symmetry operator by inspection of the reduced basis
set.

In this paper, we provide algorithms to efficiently make use of symmetry in an
operator-free manner. To do this, our methods create the subspace of
measurement basis states which would share a conserved quantity of some
commuting symmetry operator that is diagonal in that measurement basis,
without needing to explicitly create that operator. We call such a space a
symmetry-protected subspace (SPS). In the language of linear algebra, these
are invariant subspaces of the evolution operator; they are subspaces,
determined by the initial state, from which the evolution cannot escape. To
reap the benefits of symmetry we only need to find the SPS of the initial
state, not a conserved quantity, much less an explicit symmetry. However,
naively, to find an SPS we need to actually evolve the system in the full
Hilbert space, which is exponentially large in the number of qubits
(particles, spins, etc.). Sections III through V describe the formulation and
algorithms by which we avoid this exponential scaling, but here we provide a
non-technical overview.

First, note that most Hamiltonians and resulting Unitary evolution operators
are built from a number of local operators. For example, the Heisenberg-XXX
model described below consists only of nearest-neighbor interactions. So, at
some level, we have an intuition that the dynamics, thus the SPSs, should be
derivable, like the Hamiltonian itself, from a combination of local
operations, and that local operations are inherently less computationally
expensive to work with. This is indeed the case, as shown below.

Next, note that to say that a wavefunction is in a symmetry-protected or
invariant subspace is to say that it is and remains throughout dynamic
evolution a linear combination of basis states in that subspace and that
subspace alone. And if we care only about finding the subspace, we do not need
to keep track of the actual linear combination (i.e., both the basis vectors
and their amplitudes) but only the basis vectors. This “binarization” of the
evolution is critical, because it allows us to adopt a graph-theoretic
framework that is vastly more efficient for finding and searching SPSs.

This also leads to an important restriction in our work; using our methods, we
can only automate the discovery of symmetry which is diagonal in a chosen
basis. We work in the computational $Z$ basis throughout this paper, though
extensions to other bases are of course possible by rewriting the time
evolution operators in the new basis and performing the same procedure we
describe below. Our automated methods should thus be viewed as a tool to find
symmetry-protected subspaces within a given basis (if, of course, they exist),
and to potentially improve classical and noisy quantum simulations based on
those discovered subspaces. However, they still require an intelligently
guessed initial basis as a starting point. For the problems we consider in
this work, the computational basis is sufficient to derive novel results,
though more complex choices can be required in other cases.

With this caveat in mind, once the basis is chosen we create an undirected and
unweighted graph, called the state interaction graph, which describes all
possible state-to-state transitions over a single application of a unitary
evolution operator. The transitive closure of this graph fragments the Hilbert
space of the system, represented in a particular measurement basis, into a
cluster graph, whose subgraphs are each a symmetry-protected subspace.

Our main results are two classical algorithms that efficiently construct and
work within these subspaces. The first, Algorithm 1, uses “transitive closure”
on local operations to explore and explicitly construct the full SPS of an
initial state, which enables the partition of the Hilbert space into a set of
disjoint SPSs. This algorithm scales linearly in both the number of local
operators from which the global operator is constructed and the size of the
SPS, which is a huge improvement over the exponential scaling of the naive
“full evolution” approach. However, because the SPS itself can be
exponentially large (albeit with an asymptotically smaller prefactor), the
second algorithm, Algorithm 1, finds a path of local operations through a set
of SPS graphs from an initial to final (i.e., measured) state to determine if
they lie within the same SPS (and thus the final state is valid). This
algorithm scales as the number of local operations raised to a small integer
power (that can be tuned for accuracy/performance) times the length of the
path, thus completely eliminating any exponential scaling.

This paper is structured as follows. In Sec. II, we define symmetry-protected
subspaces and the three quantum systems that will serve as our “exemplars”
throughout. In Sec. III we outline our novel graph theoretical approach to
quantum simulations. Section IV discusses the algorithm for computing an
entire symmetry-protected subspace. Section V provides a more efficient
algorithm to verify if two states exist within the same symmetry-protected
subspace. Finally, we demonstrate the effectiveness of symmetry-protected
subspaces to mitigate error in quantum simulations on a classical emulator in
Sec. VI.

##  II Preliminaries

###  II.1 Symmetry-protected subspaces

Consider a quantum system undergoing unitary evolution according to the
operator $U(t)$ for a time $t$. The operator $U(t)$ can represent continuous
time evolution, but also includes other cases, such as discrete time
evolution. We say that $U(t)$ is invariant under the action of an operator $S$
if $[S,U(t)]=0$ for all times $t$. In this instance, $S$ is a symmetry
operator. For a basis in which $S$ is diagonal, states can be labeled by their
eigenvalues under the action of $S$: $S|s,b\rangle=s|s,b\rangle$, where $b$ is
some other (set of) label(s), which could represent, for example, a
computational basis state integer encoding or a many-body eigenvalue. Suppose
now that we initialize the dynamics in a state of definite $s$:
$|\psi^{s}_{0}\rangle=\sum_{b}\alpha^{s}_{b}|s,b\rangle$, and evolve under
$U(t)$. Given the commutativity of the symmetry and evolution operators, the
action of the symmetry operator on the output state is:
$S[U(t)|\psi^{s}_{0}\rangle]=U(t)S|\psi^{s}_{0}\rangle=U(t)s|\psi^{s}_{0}\rangle=s[U(t)|\psi^{s}_{0}\rangle]$,
indicating that the eigenvalue $s$ is conserved under the evolution for all
time. Our methods identify states that would share an eigenvalue under $S$
without explicitly knowing $S$. To do this, we use the notion of a symmetry-
protected subspace, also known as an invariant subspace, which we define
below.

######  Definition 1 (Symmetry-Protected Subspace).

Let $\mathcal{H}_{d}$ be a Hilbert space of dimension $d$ spanned by an
orthonormal set of basis vectors, $B(\mathcal{H}_{d})=\\{|b\rangle\\}$ . A
subspace $G\subseteq\mathcal{H}_{d}$, which is spanned by a subset of
$B(\mathcal{H}_{d})$, is a symmetry-protected subspace of unitary operator
$U(t)$ if and only if a projection onto $G$, $P_{G}=\sum_{b\in
G}\ket{b}\bra{b}$, obeys the commutation relation $[P_{G},U(t)]=0$.

Note that while our definition emphasizes the connection to symmetries, an
SPS, as defined above, is indeed an invariant subspace according to the usual
definition [19], which is simply that $\forall\ket{g}\in G,\>U(t)\ket{g}\in
G$, because if $\ket{g}\in G$,

$U(t)\ket{g}=U(t)P_{G}\ket{g}=P_{G}[U(t)\ket{g}]\in G,$ (1)

where we have used the commutativity of $P_{G}$ and $U(t)$ and the fact that
by definition the result of applying $P_{G}$ to anything is in $G$ . Another
consequence of Def. 1 is that if a particular basis state is not in $G$, then
the transition matrix element to that state from any state in $G$ is strictly
zero.

######  Lemma 0.1.

Let $\ket{g}\in G$ be an arbitrary element of a symmetry-protected (i.e.,
invariant) subspace $G$ of an evolution operator $U(t)$ and let $\ket{b}\in
B(\mathcal{H}_{d})$ be a basis vector outside of $G$, $\ket{b}\notin G$. Then,
$\bra{b}U(t)\ket{g}=0$ for any time $t$.

###### Proof.

Using Def. 1,
$\bra{b}U(t)\ket{g}=\bra{b}U(t)P_{G}\ket{g}=\bra{b}P_{G}U(t)\ket{g}=[P_{G}\ket{b}]^{\dagger}U(t)\ket{g}=0$,
since the projection operator onto $G$ annihilates states outside of $G$. ∎

In Sec. IV we present an algorithm that prescriptively constructs a subspace,
denoted $G_{\ket{\psi_{0}}}$, of a particular initial state, $\ket{\psi_{0}}$.
We also show in the corresponding theorem, Thm. 1, that subspaces so
constructed share the property described in Lem. 0.1, which relies on Def. 1
for its proof. Hence, $G_{\ket{\psi_{0}}}$ constructed according to the
procedure in Sec. IV are symmetry-protected subspaces.

Similar to that of the symmetry operator $S$, the commutation relation
involving $P_{G}$ also identifies a dynamical invariance, since for any
$|g\rangle\in G$, $P_{G}[U(t)|g\rangle]=[U(t)|g\rangle]$. That is, the state
$U(t)|g\rangle$ is an eigenvector of the projector with eigenvalue $1$.
Therefore, while the projector $P_{G}$ does not identify conserved symmetry
eigenvalues, it does indicate when such an eigenvalue exists. In Sec. III, we
will show how discovering $P_{G}$, rather than $S$ directly, empowers our
algorithms to discover underlying invariances.

###  II.2 Quantum simulations

The general form for unitary evolution operators we assume is

$U(t)=\mathcal{O}_{op}\Bigg{[}\prod_{j=1}^{p}\prod_{i=1}^{m}U_{i}(\tau_{j})\Bigg{]},$
(2)

where $t$ is the total duration of evolution, $\mathcal{O}_{op}$ denotes some
operator ordering, such as time ordering, and $p$ is the number of time-steps
used to evolve to $t$. Each $U_{i}(\tau_{j})$ acts locally on $k$-qubits and
is parameterized by a real time coordinate $\tau_{j}$. Evolution unitaries of
the form in Eq. (2) can evolve over discrete time, where each $\tau_{j}$ is
finite, or discretized-continuous time, where each $\tau_{j}$ is (ideally
infinitessimally) small in order to minimize Trotter error.

In the instance where a $k$-local Hamiltonian, $H=\sum_{i}h_{i}$, is known,
Eq. (2) results from the local dynamics governed by the $h_{i}$ via the
Trotter-Suzuki formula [20], and $p$ corresponds to the number of Trotter
steps. We refer to the operator $U$ as the relevant quantum system and examine
the dynamics and associated symmetry-protected subspaces of three exemplary
systems.

####  II.2.1 Heisenberg-XXX

The one-dimensional Heisenberg-XXX model for $n$ spin-$1/2$ particles with
nearest-neighbor interactions is given by the following Hamiltonian:

$H=\sum_{i=0}^{n-2}X_{i}X_{i+1}+Y_{i}Y_{i+1}+Z_{i}Z_{i+1},$ (3)

where $X_{i},Y_{i},$ and $Z_{i}$ are Pauli operator acting on spin $i$. The
model conserves total spin in the $Z$ basis, represented by the operator

$S^{z}=\sum_{i=0}^{n-1}Z_{i},$ (4)

as well as the correspondingly-defined operators $S^{x}$ and $S^{y}$. Quantum
simulation of the Heisenberg model on a digital quantum processor can be
achieved via exponentiation and Trotterization of Eq. (3). In such quantum
simulation experiments, one usually picks one or a number of qubit bases in
which to measure. The symmetry operators $S^{x,y,z}$ can be used to mitigate
errors in the all-qubit $X,Y,Z$ measurement bases, respectively, via post-
selection. In the context of classical simulation of the XXX model, the
symmetries can be used to constrain the number of basis states included in the
dynamics.

####  II.2.2 $T_{6}$ quantum cellular automata

The one-dimensional $T_{6}$ quantum cellular automata (QCA) rule has recently
come to interest within the context of quantum complexity science as a
dynamical small-world mutual information network generator [21] and a QCA
Goldilocks rule [22]. Its discrete-time unitary update can be derived from a
parent Hamiltonian, but it is more natural to define the system by specifying
the simulation unitary for a discrete time $t=p$, directly:

$\begin{split}U(T_{6};t)&=\mathcal{O}_{op}^{time}\Bigg{[}\prod_{j=1}^{p}\prod_{i=3,5,\ldots}^{n-1}U_{i}(\tau_{j})\prod_{i=2,4,\ldots}^{n-1}U_{i}(\tau_{j})\Bigg{]}\\\
U_{i}(\tau_{j})&=\sum_{\alpha,\beta=0}^{1}P_{i-1}^{(\alpha)}\otimes(H_{i})^{\delta_{\alpha+\beta,1}}\otimes
P_{i+1}^{(\beta)},\end{split}$ (5)

where $P^{(\alpha)}_{i}=\ket{\alpha_{i}}\bra{\alpha_{i}}$ for $\alpha=0,1$ is
the projection operator onto the corresponding state of qubit $i$, $H_{i}$ is
the Hadamard operator, and $\delta_{\alpha+\beta,1}$ is the Kronecker delta
function. At each time step, a Hadamard is applied to a qubit only if exactly
one of its neighbors is in the $\ket{1}$ state (i.e., $\alpha+\beta=1$) and
does nothing otherwise. It has a known $Z$ basis symmetry related to domain-
wall conservation:

$S=\sum_{i=0}^{n}Z_{i}Z_{i+1},$ (6)

where it should be understood that indices $i=1,\ldots,n$ refer to dynamical,
computational qubits while the indices $i=0$ and $i=n+1$ refer to non-
dynamical qubits fixed to the $\ket{0}$ state.

####  II.2.3 $F_{4}$ quantum cellular automata

The one-dimensional $F_{4}$ QCA with nearest- and next-nearest-neighbor
connectivity is another Goldilocks rule [22]. It is also most easily specified
by its simulation unitary for discrete-time duration $t=p$:

$\begin{split}U(F_{4};t)=\mathcal{O}_{op}^{time}\Bigg{[}\prod_{j=1}^{p}U(\tau_{j})\Bigg{]},\end{split}$
(7)

where if the time step index, $j$, is even

$U(\tau_{j=\text{even}})=\prod_{i=2,5,8,\ldots}U_{i}(\tau_{j})\prod_{i=3,6,9,\ldots}U_{i}(\tau_{j})\prod_{i=4,7,10,\ldots}U_{i}(\tau_{j})$
(8)

and if $j$ is odd, then

$U(\tau_{j=\text{odd}})=\prod_{i=3,6,9,\ldots}U_{i}(\tau_{j})\prod_{i=2,5,8\ldots}U_{i}(\tau_{j})\prod_{i=4,7,10,\ldots}U_{i}(\tau_{j}).$
(9)

In either case,

$U_{i}(\tau_{j})=\sum_{\alpha,\beta,\gamma,\omega=0}^{1}P_{i-2}^{\alpha}P_{i-1}^{\beta}(H_{i})^{\delta_{\alpha+\beta+\gamma+\omega,2}}P_{i+1}^{\gamma}P_{i+2}^{\omega}.$
(10)

Equation (10) applies a Hadamard to a qubit if exactly two out of its
neighbors or next-nearest neighbors are in the $\ket{1}$ state. There are no
analytically known symmetries for this rule. As shown in Fig. 3 in Sec. IV.3,
our methods discover previously-unknown, symmetry-protected subspaces,
indicating a hitherto hidden symmetry of the system.

FIG. 1: Example construction of the state interaction graph
$D_{U_{\text{Trot}}}$ and corresponding symmetry-protected subspaces for the
one-dimensional, four-site hopping unitary in Eq. (12). Multi-qubit states are
ordered as $\ket{q_{0}q_{1}q_{2}q_{3}}$. One Trotter approximation for the
unitary, $U_{\text{Trot}}$, is given in Eq. (13), where it should be
understood that right-most operators act first. For graphical clarity we omit
loops, with the understanding they are always implied. (a) State interaction
graph $D_{U_{\text{Trot}}}$. Nodes represent the different 4-bit strings, and
edges occur where matrix elements of $U_{\text{Trot}}$ are nonzero. Note the
presence of five disconnected subgraphs. In practice, we treat each of these
edges as undirected. (b) Transitive closure of $D_{U_{\text{Trot}}}$ resulting
in the closed state interaction graph, which in this case corresponds exactly
to the state interaction graph for the hopping unitary:
$D_{U_{\text{Hop}}}=D^{+}_{U_{\text{Trot}}}$. Each disconnected subgraph in
$D_{U_{\text{Trot}}}$ has become a disconnected complete graph in
$D^{+}_{U_{\text{Trot}}}$. (c) Each complete graph in
$D_{U_{\text{Hop}}}=D^{+}_{U_{\text{Trot}}}$ corresponds to a symmetry-
protected subspace $G_{s}$. In this example, the subspace indices correspond
to eigenvalues under particle number, $S=\sum_{i=0}^{3}a^{\dagger}_{i}a_{i}$,
conservation.

##  III Graph Theory Approach to Quantum Simulations

In this section, we show how graph theory coupled with transitive closure
discovers symmetry-protected subspaces. We describe how to “binarize”
interactions between basis states through the dynamical system by discarding
amplitude and phases to solely highlight if such an interaction exists,
describing these interactions as transitive relations, and finding the states
connected by a transitive relation to form the SPSs of the system.

Directly, this method still requires unitary matrix multiplication in Hilbert
space to establish basis state interactions through the dynamical system.
Thus, to apply our methods, we need an efficient way to describe state-to-
state interactions. We do this by creating a structure we call a string edit
map $\mathcal{L}_{U}$ for a unitary operator $U$, which relies on the
observation that quantum systems are typically structured by local
interactions. This string edit map returns the basis states available through
any operations of a single local unitary operator on a single basis state in
near-constant time complexity, allowing us to inexpensively find basis vectors
available in local quantum dynamics.

This section will proceed as follows: in Sec. III.1, we first define the
concept of a state interaction graph, whose edges indicate nonzero amplitudes
on transitions between measurement basis states in the quantum simulation.
Next, in Sec. III.2 we show how transitive closure on this graph creates a
cluster graph, which we call the closed state interaction graph, whose
complete subgraphs are symmetry-protected subspaces. Finally, in Sec. III.3 we
show how this closed interaction graph can efficiently return the set of all
measurement basis states seen by a unitary operator on a single state via the
construction of the string edit map $\mathcal{L}_{U}$ for a unitary operator
$U$.

###  III.1 State interaction graph

We begin by defining the state interaction graph $D_{U}$ for a unitary
operator $U$ over a Hilbert space basis $B(\mathcal{H}_{d})$ and showing how
to construct it. For now we leave the form of $U$ general and will specify
particular forms when necessary.

######  Definition 2 (State Interaction Graph).

Given a basis $B$ and a unitary operator $U$, define a vertex set $V\equiv B$
and an undirected edge set
$E\equiv\\{(\ket{b}\leftrightarrow\ket{b^{\prime}})\>\>\forall\>\>\ket{b},\ket{b^{\prime}}\in
B\>\>|\>\>\bra{b^{\prime}}U\ket{b}\neq 0\\}$. Then, the state interaction
graph is defined by the ordered tuple $D_{U}\equiv(V,E)$. In other words, the
basis states of $B$ are assigned to vertices (nodes) in the interaction graph
and edges are created between vertex states only where the matrix element of
the evolution unitary between the two states is nonzero.

Strictly speaking, $D_{U}$ should be a directed graph, where an edge points
from $\ket{b}$ to $\ket{b^{\prime}}$ if $\bra{b^{\prime}}U\ket{b}\neq 0$ and
from $\ket{b^{\prime}}$ to $\ket{b}$ if $\bra{b}U\ket{b^{\prime}}\neq 0$.
However, for a symmetry operator $S$ and time-dependent simulation unitary
$U(t)$, one can show that $[S,U(t)]=0\iff[S,U^{\dagger}(t)]=0$, meaning that
symmetries of the evolution operator, and their associated protected
subspaces, are invariant under time-reversal. To remain consistent with this
observation, we treat every directed relationship
$\bra{b^{\prime}}U\ket{b}\neq
0\Rightarrow(\ket{b}\rightarrow\ket{b^{\prime}})$ as an undirected edge
$(\ket{b}\leftrightarrow\ket{b^{\prime}})$. This treatment is equivalent to
assuming that the true, directed state interaction graph corresponding to Def.
2 always has a cycle that leads back to every node, such as is the case in
Fig. 1(a). Formally, we assume that if an edge
$(\ket{b}\rightarrow\ket{b^{\prime}})\in D_{U}$, there exists a path
$\\{\ket{b^{\prime}}\rightarrow\dots\rightarrow\ket{b}\\}\subseteq D_{U}$.
This global cyclicity assumption enables us to use the notion of transitive
closure in Sec. III.2, and subsequently, in an uncomplicated manner that
respects the time-reversal invariance of the resulting subspaces. It is
justified on two points: 1.) Most simulation unitaries have a very regular,
repetitive structure so that directed, acyclic state interaction graphs are
likely only to arise in extremely pathological instances, and 2.) The failure
of the global cyclicity assumption will only ever result in the artificial
enlargement of a symmetry-protected subspace, and while such a failure leads
to underconstrained subspaces, which is bad for the efficacy of e.g., post-
selection, it will never result in the corruption of simulation fidelity by
overconstraining or throwing out good simulation data.

To construct $D_{U}$, we use the following steps. First, choose a set of
Hilbert space basis vectors $B(\mathcal{H}_{d})=\\{\ket{b}\\}$. Any symmetry-
protected subspaces that are discovered must be formed by the basis vectors of
this basis. In the context of quantum simulation, $B$ dictates the basis in
which a quantum computer will be measured. For example, parallel readout in
the computational $Z$ basis will result in bit strings,
$B_{\text{comp.}}=\\{\ket{0\ldots 00},\ket{0\ldots 01},\ldots,\ket{1\ldots
11}\\}$, which are Pauli $Z$-string eigenvectors. In the context of classical
simulation, the basis furnishes a representation for the $d$-dimensional
vector of complex amplitudes that stores the evolving many-body wavefunction.
When the evolution operator is applied to a basis vector for a single time-
step the resulting state $\ket{\psi(\tau)}$ has a basis vector decomposition

$\ket{\psi(\tau)}=U(\tau)\ket{b}=\sum_{b^{\prime}\in
B}\ket{b^{\prime}}\bra{b^{\prime}}U(\tau)\ket{b}=\sum_{b^{\prime}\in
B}\alpha_{b^{\prime}}(\tau)\ket{b^{\prime}},$ (11)

where we suppress the time-ordering subscript $j$ in $\tau_{j}$ for
simplicity. With this decomposition for each $\ket{b}\in B$, we create the
state interaction graph $D_{U}$ of the operator $U(\tau)$: for each pair
$\ket{b},\ket{b^{\prime}}$ such that $\alpha_{b^{\prime}}(\tau)\neq 0$ in Eq.
(11), we add an edge $(\ket{b}\leftrightarrow\ket{b^{\prime}})$ to $D_{U}$. If
one has a direct $d\times d$ matrix representation of $U(\tau)$ on-hand, then
the adjacency matrix for $D_{U}$ can be read off directly as
$\mathcal{A}_{U}=\text{bit}[U(\tau)]+\text{bit}[U^{\dagger}(\tau)]$, where if
an entry in $U(\tau)$ or $U^{\dagger}(\tau)$ becomes non-zero for any value of
$t$, its complex value is replaced by $1$ under the operation
$\text{bit}[\ldots]$. This adjacency matrix can be formed, for example, by
directly exponentiating a $d\times d$ Hamiltonian matrix with a time
parameter. In practice, however, constructing or storing an entire evolution
unitary in memory is costly, since the size of Hilbert space grows
exponentially in the number of qubits, or spins, in the simulation: $d=2^{n}$.
Indeed, one of the main advantages of digital quantum simulation is the
ability to break global evolution unitaries into sequences of local unitaries,
at the expense of introducing error, which are then implemented as quantum
gates. Therefore, being able to extract symmetry-protected subspaces from
consideration of local operations, rather than from the global unitary they
may approximate, is of clear benefit.

Towards this end, Fig. 1(a) shows the state interaction graph
$D_{U_{\text{Trot}}}$ for one, potentially very bad depending on $\theta$,
Trotter approximation to the one-dimensional, four-qubit hopping unitary

$U_{\text{Hop}}(\theta)=e^{i\theta\sum_{i=0}^{2}(X_{i}X_{i+1}+Y_{i}Y_{i+1})/2}.$
(12)

We take the Trotterization to be

$U_{\text{Trot}}(\theta)=\text{iSWAP}_{01}(\theta)\times\text{iSWAP}_{12}(\theta)\times\text{iSWAP}_{23}(\theta),$
(13)

where it should be understood that rightmost operators are applied first, and
multi-qubit states are ordered as $\ket{q_{0}q_{1}q_{2}q_{3}}$. Notice that
each $\text{iSWAP}_{i,i+1}(\theta)$ is parameterized by an arbitrary $\theta$,
and as such we expect each iSWAP operation to be a fractional operation that
leaves some residual state behind, i.e., the operator has an identity
component. The system is small enough that the state interaction graph can be
checked by hand in this case, and the main observation to be made is that it
is comprised of four, disjoint subgraphs, each only containing transitions
between states of fixed particle number (i.e., number of $\ket{1}$s), and that
all nodes in each subgraph have a path to all other nodes in the subgraph. The
zero- and four-particle states are isolated, while the one-, two-, and three-
particle states form directed, incomplete, isolated subgraphs. As we will see,
this “incompleteness” feature is a pathology of the Trotter approximation
which will be rectified in Sec. III.2 via transitive closure. It is also worth
noting that constructing the state interaction graph using the Trotter-
approximated unitary is not yet useful, since on a classical computer it
currently still requires the storage and evolution of a $2^{n}$-dimensional
wavefunction. We will demonstrate the utility of constructing state
interaction graphs from component $k-$local operators in Sec. III.3.

###  III.2 Defining symmetry-protected subspaces with transitive closure

Suppose we have states $\ket{b}$ and $\ket{b^{\prime\prime}}$, such that
$\bra{b^{\prime\prime}}U\ket{b}=0$ and $\bra{b^{\prime\prime}}UU\ket{b}\neq
0$. This requires a transitive relation:
$\bra{b^{\prime\prime}}UU\ket{b}=\sum_{b^{\prime}\in
B}\bra{b^{\prime\prime}}U\ket{b^{\prime}}\bra{b^{\prime}}U\ket{b}$, because
$\ket{b}$ must first transition to an intermediate state $\ket{b^{\prime}}$ to
reach its final destination at $\ket{b^{\prime\prime}}$. Therefore, in our
state interaction graph $D_{U}$, as defined in Def. 2, there are edges
$(\ket{b}\rightarrow\ket{b^{\prime}}),(\ket{b^{\prime}}\rightarrow\ket{b^{\prime\prime}})\in
D_{U}$, which will have the same transitive relation encoded in the path
$\\{\ket{b}\rightarrow\ket{b^{\prime}}\rightarrow\ket{b^{\prime\prime}}\\}\subseteq
D_{U}$. We use this duality to make the assumption that if the edges
$(\ket{b}\rightarrow\ket{b^{\prime}}),(\ket{b^{\prime}}\rightarrow\ket{b^{\prime\prime}})\in
D_{U}$, the transition amplitude $\bra{b^{\prime\prime}}UU\ket{b}\neq 0$.

There are cases where the amplitudes of the states cancel, due to destructive
interference, and break this transitive property on the level of individual
basis state to basis state interactions. By ignoring the amplitudes of the
basis states, we run the risk of including states in the subspace which would
be removed via destructive interference. This risk comes with the benefit of
efficiently knowing which states are reachable in the quantum simulation, and
for our applications it does not add any error to a simulation, as it does not
break the underlying commuting subspace projection operator $P_{G}$; the
subspaces are simply not as restrictive as they could be. See Appendix A for a
proof. This also allows one to define the symmetry-protected subspaces to be
for any parameterization of the simulation unitary.

The transitive property exists for every state $\ket{b}$ in $D_{U}$, so we can
take the transitive closure of the state interaction graph to create a closed
state interaction graph $D^{+}_{U}$. The transitive closure of an edgeset $E$
is a transitively closed edgeset $E^{+}$, where every pair of states
$\ket{b},\ket{b^{\prime}}\in E$ which can be associated by any transitive
relation, in other words can be connected by a path
$\\{\ket{b}\rightarrow\dots\rightarrow\ket{b^{\prime}}\\}\subseteq E$, has an
edge $(\ket{b}\rightarrow\ket{b^{\prime}})\in E^{+}$ [23].

######  Definition 3 (Transitively Closed State Interaction Graph).

Let $D_{U}=(V,E)$ be a potentially non-closed state interaction graph for
unitary $U$ and basis $B$. Define $V^{+}\equiv V$ to be the closed state
interaction graph vertex (node) set and $E^{+}$ to be the state interaction
graph edge set. An edge, $(\ket{b}\leftrightarrow\ket{b^{\prime}})\in E^{+}$,
exists in this edge set if and only if there is a path between $\ket{b}$ and
$\ket{b^{\prime}}$ in $D_{U}$,
$\\{\ket{b}\leftrightarrow\dots\leftrightarrow\ket{b^{\prime}}\\}\subseteq E$.
The transitively closed state interaction graph is then defined as
$D_{U}^{+}\equiv(V^{+},E^{+})$.

Because the original interaction graph represents single-operator state-to-
state transitions, any two basis states which can discover each other through
the quantum evolution have an edge in $D^{+}_{U}$; in other words,

$\exists\>t\;s.t.\;\bra{b^{\prime}}U(t)\ket{b}\neq
0\Rightarrow(\ket{b}\leftrightarrow\ket{b^{\prime}})\in D^{+}_{U},$ (14)

for some time/operator exponent $t$. The transitive closure of an undirected,
unweighted graph is a cluster graph, or a set of complete subgraphs; as
discussed at the end of Sec. III.1, we treat the graph $D_{U}$ as undirected
just for this purpose. The transitive closure of the state interaction graph
for the Trotterized unitary in Fig. 1(a) can be seen in Fig. 1(b), which in
this case turns out to be state interaction graph for the original hopping
unitary, or $D^{+}_{U_{\text{Trot}}}=D_{U_{\text{Hop}}}$. As with
$D_{U_{\text{Trot}}}$, within each complete subgraph, total particle number is
conserved.

If the edge $(\ket{b}\leftrightarrow\ket{b^{\prime}})\in D^{+}_{U}$, then
$\ket{b}$ and $\ket{b^{\prime}}$ share a conserved quantity of the underlying
unknown symmetry. This cluster graph structure also makes it apparent that if
a wavefunction is initialized as a linear combination of vectors in one
subgraph $G$ of $D^{+}_{U}$: $\ket{\psi_{0}}=\sum_{b\in G}\alpha_{b}\ket{b}$,
it will remain in that subgraph:

$\ket{\psi(t)}=U(t)\ket{\psi_{0}}=\sum_{b^{\prime}\in
G}\alpha_{b^{\prime}}(t)\ket{b^{\prime}}\;\forall\>t.$ (15)

Therefore, each complete subgraph $G$ within $D^{+}_{U}$ is a symmetry-
protected subspace. We will formally prove that transitive closure on the
state interaction graph can give an SPS with Thm. 1 in Sec. IV.1.

$D^{+}_{U}$ can be represented as a list of disjoint sets of nodes, where each
set has implied all-to-all connectivity. This set construction can be seen for
the hopping unitary in Fig. 1(c). Here, the index, $s$, of each subset,
$G_{s}$ counts the number of conserved particles. Formally, we can define the
associated symmetry-protected subspaces by constructing their projection
operators according to Def. 1: $P_{G}=\sum_{b\in G}\ket{b}\bra{b}$.

###  III.3 Basis state string edit map

We now have a method to identify symmetry-protected subspaces using the
language of graph theory. However, actually computing these subspaces still
requires the construction and manipulation of vectors and matrices in an
exponentially large Hilbert space. Recall, though, that the systems of
interest are defined by Hamiltonians composed of local operations
$U=\prod_{i}U_{i}$ where each $U_{i}$ is k-local, meaning it only involves $k$
of the $n$ total qubits in the system, and where in general we will have $k\ll
n$. In this section we describe our mechanism for using this fact to build up
SPSs efficiently with what we call the basis string edit map. This map enables
computation of subspaces using only $k$-local operations on basis state
vectors, so the computational complexity of operation with this map scales
with $k$ instead of $n$.

######  Definition 4 (Basis String Edit Map).

Let $U_{i}(t)$ be a unitary operator that acts for a time $t$ non-trivially on
$k$ of $n$ qubits, $Q_{k}(i)\subseteq\\{q_{0},\ldots,q_{n-1}\\}$. The basis
string edit map $\mathcal{L}_{U_{i}}$ maps a basis state $\ket{b}$ to the set
of basis states $\\{\ket{b^{\prime}}\\}$ to which $\ket{b}$ can evolve after
an arbitrary amount of time under $U_{i}(t)$. Formally,
$\mathcal{L}_{U_{i}}(\ket{b})=\\{\ket{b^{\prime}}\>|\>\exists\,t\,\bra{b^{\prime}}U_{i}(t)\ket{b}\neq
0\\}$.

We can apply this construction to any unitary operator, including any
$k$-local Trotter decomposition. Given a unitary operator which is a product
of local operators $U=\prod_{i}U_{i}$, we form the set of local operators used
in the Trotter decomposition, $\\{U_{i}\\}$. Then a basis string edit map can
be formed for any subset of operators from $\\{U_{i}\\}$, as long as every
operator $U_{i}\in U$ is included in at least one string edit map. We will use
the decomposition

$\mathbb{L}\equiv\\{\mathcal{L}_{U_{i}}:U_{i}\in U\\},$ (16)

which has one string edit map for each local operator in the Trotter
decomposition.

In Def. 4 we have deliberately left out the exact space upon which
$\mathcal{L}_{U_{i}}$ acts. When operating on states in a basis $B$ with
$\mathcal{L}_{U_{i}}$, when $\text{dim}(\mathcal{L}_{U_{i}})<\text{dim}(B)$,
we will call $\mathcal{L}_{U_{i}}$ a “substring edit map”, and when
$\text{dim}(\mathcal{L}_{U_{i}})=\text{dim}(B)$, we refer to it as a “string
edit map” or “full string edit map”.

For example, in the hopping unitary, $U_{\text{Hop}}$ given in Eq. (12), we
can define $\mathcal{L}_{U_{\text{Hop}}}$ which would give
$\mathcal{L}_{U_{\text{Hop}}}(\ket{0111})=G_{3}$ where $G_{3}$ is in Fig.
1(c). More usefully, we can take the Trotterization of Eq. (12) given by Eq.
(13), and define string edit maps for the local iSWAP operators,
$\mathbb{L}_{\text{Trot}}=\\{\mathcal{L}_{\text{iSWAP}_{0,1}(\theta)},\mathcal{L}_{\text{iSWAP}_{1,2}(\theta)},\mathcal{L}_{\text{iSWAP}_{2,3}(\theta)}\\}$.
In Eq. (13), each local operator $U_{i}=\text{iSWAP}_{i,i+1}(\theta)$ would
generate the corresponding edit map
$\mathcal{L}_{\text{iSWAP}_{i,i+1}(\theta)}$, which operate as, e.g.,
$\mathcal{L}_{\text{iSWAP}_{1,2}(\theta)}(\ket{0\textbf{10}0})=\\{\ket{0\textbf{01}0},\ket{0\textbf{10}0}\\}$.
Here we have written the operators as acting on the full $2^{n}$-dimensional
Hilbert space and highlighted in bold the qubits that are part of the $2^{k}$
dimensional subset of this space upon which
$\mathcal{L}_{\text{iSWAP}_{1,2}(\theta)}$ acts. This k-local string edit map
does not require information about any states besides those at the relevant
indices, $1$ and $2$ in this case.

Algorithm 1 creates $\mathcal{L}_{U_{i}}$ by taking transitive closure of the
adjacency matrix $\mathcal{A}_{U_{i}}$ of $D_{U_{i}}$ via boolean matrix
multiplication [24]. This algorithm requires $O(2^{3k})$ time to compute
$\mathcal{L}_{U_{i}}$, $O(2^{k})$ space to store it, and $O(1)$ time to use
$\mathcal{L}_{U_{i}}$.

Algorithm 1 Create a string edit map $\mathcal{L}_{U_{i}}$

Local unitary $U_{i}$, orthonormal basis $B_{i}=\\{\ket{b}\\}$

$\text{Dim}(B_{i})=\text{Dim}(U_{i})=2^{k}$ and
$\text{Col}(U_{i})=\text{Col}(B_{i})$

$\mathcal{A}_{U_{i}}\leftarrow\text{bit}[U_{i}]+\text{bit}[U_{i}^{\dagger}]$
$\triangleright$ Add $U_{i}^{\dagger}$ to make the adjacency matrix undirected

$\mathcal{A}_{U_{i}}^{\prime}\leftarrow\text{bit}[\mathcal{A}_{U_{i}}^{2}]$

while $\mathcal{A}_{U_{i}}^{\prime}\neq\mathcal{A}_{U_{i}}$ do

$\mathcal{A}_{U_{i}}\leftarrow\mathcal{A}_{U_{i}}^{\prime}$

$\mathcal{A}_{U_{i}}^{\prime}\leftarrow\text{bit}[\mathcal{A}_{U_{i}}^{2}]$

end while

for $\ket{b}\in B_{i}$ do

for $\ket{b^{\prime}}\in B_{i}$ do

if $\mathcal{A}_{U_{i}}[\ket{b},\ket{b^{\prime}}]=1$ then

$\mathcal{L}_{U_{i}}(\ket{b})\leftarrow\mathcal{L}_{U_{i}}(\ket{b})\cup\ket{b^{\prime}}$

end if

end for

end for

return $\mathcal{L}_{U_{i}}$

Algorithm 1 is computationally trivial to compute for small unitary operators
($k=2$ in the iSWAP example), but very expensive for the large unitary
operators encountered in quantum simulations. Throughout the rest of this
paper, we will use $\mathcal{L}_{U_{i}}$ defined on small $k$ to compute our
subspaces in order to keep a small overhead (thus, in our terminology, we will
always be talking about “substring edit maps”).

The object $\mathbb{L}$ lets us create a version of the state interaction
graph from Def. 2 which we call the string interaction graph $D_{\mathbb{L}}$;
this construction compactly shows the state-to-state interactions in the
algorithms presented in Sec. IV and V. The mapping itself is defined as
$\mathbb{L}(\ket{b})\equiv\\{\mathcal{L}_{i}(\ket{b})\>:\>\mathcal{L}_{i}\in\mathbb{L}\\}$.

######  Definition 5 (String Interaction Graph).

Let $D_{\mathbb{L}}$ be a graph defined by an ordered tuple
$D_{\mathbb{L}}\equiv(V_{\mathbb{L}},E_{\mathbb{L}})$ and
$\mathbb{L}\equiv\\{\mathcal{L}_{U_{i}}\>:\>U_{i}\in U\\}$ for some
Trotterized unitary $U=\prod U_{i}$. The vertex set is then given by
$V_{\mathbb{L}}\equiv B(\mathcal{H}_{d})$ and the edgeset is given by
$E_{\mathbb{L}}\equiv\\{(\ket{b}\leftrightarrow\ket{b^{\prime}})\>:\>\ket{b},\ket{b^{\prime}}\in
B\>\mathrm{and}\>\ket{b^{\prime}}\in\mathbb{L}(\ket{b})\\}$.

Notice that, while edges in $D_{U}$ can capture the action of multiple
operators at once, each edge in $D_{\mathbb{L}}$ is the action of only a
single operator; therefore, if $D_{U}\equiv(V,E)$ is the state interaction
graph and $D_{\mathbb{L}}\equiv(V,E_{\mathbb{L}})$ is the string interaction
graph, then $E_{\mathbb{L}}\subseteq E$. Later, we will show with the proof in
Appendix B for Thm. 1 that despite this inequality, the transitive closure of
the graph $D_{\mathbb{L}}^{+}$, is equivalent to the transitively closed state
interaction graph $D_{\mathbb{L}}^{+}\equiv D_{U}^{+}$.

FIG. 2:  This figure shows how we can incrementally build the symmetry-
protected subspace for $U=\prod_{i=0}^{2}\text{iSWAP}_{i,i+1}(\theta)$. (a) A
single un-closed symmetry-protected subspace, or a subgraph of
$D_{\mathbb{L}}$ for an iSWAP network. All operations
$\mathbb{L}_{\text{Trot}}=\\{\mathcal{L}_{\text{iSWAP}_{0,1}(\theta)},\mathcal{L}_{\text{iSWAP}_{1,2}(\theta)},\mathcal{L}_{\text{iSWAP}_{2,3}(\theta)}\\}$
which are non-identity on each node are shown as an edge. Notice the
similarities between this graph and the $s=2$ subgraph of Fig. 1(a): both have
the same vertex set, but the edgeset of $D_{\mathbb{L}}$ is a subset of the
edgeset in $D_{U}$. (b) Follow the recursion relation in Eq. (17) to
iteratively build the symmetry-protected subspace $G_{\ket{1100}}$, starting
from $\ket{\psi_{0}}=\ket{1100}$. Even though the graph in (a) is not
equivalent to the graph in Fig. 1(a), their transitively closed graphs are
equivalent, as can be seen by the vertices covered by the red line.

##  IV Algorithm: creating the symmetry-protected Subspace of an Initial State

We now turn to the construction of the entire symmetry-protected subspace of a
given initial state using $k$-local substring edit maps and transitive
closure. The algorithm by which we do so works by establishing a recurrence
relation for computing new states, associated by symmetry-protection to the
initial state, in the simulation. This recurrence relation furnishes an
efficient way to compute the transitive closure of select subgraphs of the
entire Hilbert space, with no extraneous information. For a unitary simulation
operator decomposed into enumerated local operations,
$U(\tau_{j})=\prod_{i=1}^{m}U_{i}(\tau_{j})$, where each $U_{i}(\tau_{j})$ is
$k$-local, we will show that a symmetry-protected subspace
$G_{\ket{\psi_{0}}}$ of the initial state $\ket{\psi_{0}}$ costs
$O((m+1)\times|G_{\ket{\psi_{0}}}|)$ to compute and $O(|G_{\ket{\psi_{0}}}|)$
to store with a breadth-first search [25]. We will give the algorithm for the
case where the initial state $\ket{\psi_{0}}$ is a single measurement basis
state (i.e., a product state). If $\ket{\psi_{0}}$ is a linear combination of
measurement basis vectors, the algorithm can be repeated for each basis vector
in the sum; this does not impact the asymptotic performance of the algorithm,
as it only makes the computed subspace bigger.

We enumerate every state in a symmetry-protected subspace by transitively
closing subgraphs created by the local basis substring edit maps established
in Sec. III.3, which return a set of basis states evolved to by their
corresponding unitary operators in $O(1)$ when operating on a single basis
state. We recursively build the subspace by checking the set of substring edit
maps $\mathbb{L}\equiv\\{\mathcal{L}_{U_{i}}\>:\>U_{i}\in U\\}$ on each new
state, until none are added. This process can be seen as the transitive
closure of a subgraph of the graph $D_{\mathbb{L}}$. For $U_{\text{Trot}}$ in
Eq. (13), Fig. 2(a) shows an example of the subgraph of $D_{\mathbb{L}}$
corresponding to the action of each
$\mathcal{L}_{\text{iSWAP}_{i,i+1}(\theta)}\in\mathbb{L}_{\text{Trot}}$
starting from the initial state $\ket{\psi_{0}}=\ket{1100}$ (self-edges are
ignored as elsewhere in the manuscript).

###  IV.1 Recurrence relation

To find the symmetry-protected subspace of an initial state, we begin by
computing the string edit map for each $k$-local unitary,
$\mathbb{L}\equiv\\{\mathcal{L}_{U_{i}}\>:\>U_{i}\in U\\}$. Next, we check the
set of measurement basis strings generated by operating with each substring
edit map on the initial state, notated
$\mathbb{L}(\ket{\psi_{0}})\equiv\\{\mathcal{L}_{i}(\ket{\psi_{0}})\>:\>\mathcal{L}_{i}\in\mathbb{L}\\}$.
We define the set
$T_{\ket{\psi_{0}}}^{1}\equiv\\{\ket{\psi_{0}}\\}\>\cup\>\mathbb{L}(\ket{\psi_{0}})$.
Then, for each new state $\ket{\phi}\in\mathbb{L}(\ket{\psi_{0}})$, we check
operations under the substring edit maps:
$T_{\ket{\psi_{0}}}^{2}=\mathbb{L}(T_{\ket{\psi_{0}}}^{1})\equiv\\{\mathbb{L}(\ket{\phi})\>:\>\ket{\phi}\in\mathbb{L}(\ket{\psi_{0}})\\}$.
This process repeats until no new states are found through the following
recurrence relation.

$\begin{cases}T_{\ket{\psi_{0}}}^{1}\leftarrow\\{\ket{\psi_{0}}\\}\cup\mathbb{L}(\ket{\psi_{0}})&\text{base
case}\\\ T_{\ket{\psi_{0}}}^{i+1}\leftarrow
T_{\ket{\psi_{0}}}^{i}\cup\mathbb{L}(T_{\ket{\psi_{0}}}^{i})&\text{recursive
case}\\\ T_{\ket{\psi_{0}}}^{i+1}=T_{\ket{\psi_{0}}}^{i}&\text{stop
condition}\end{cases}$ (17)

The stop condition activates if no new states are found, that is, when
additional operations drawn from $\mathbb{L}$ do not unveil any new states.
Steps 0-4 in Fig. 2(b) show how the recurrence relation manifests for the
input state $\ket{\psi_{0}}=\ket{1100}$ and the set of substring edit maps
generated from the Trotterization in Eq. (13). We then define the symmetry-
protected subspace $G_{\ket{\psi_{0}}}$ to which the state $\ket{\psi_{0}}$
belongs via

$G_{\ket{\psi_{0}}}\equiv T_{\ket{\psi_{0}}}^{i+1},$ (18)

where the definition “$\equiv$” in Eq. (18) should be taken to mean “all basis
states in $T_{\ket{\psi_{0}}}^{i+1}$ viewed as nodes in a complete graph”. We
can create the Kleene Closure of the set of substring edit maps, denoted
$\mathbb{L}^{\star}$, which is the set of all finite concatenations of
substring edit maps, including the identity. Any arbitrary string of substring
edit maps $\mathbb{L}^{\star}$ applied to $\ket{\psi_{0}}$ will result in a
state in $G_{\ket{\psi_{0}}}$ by its definition. Hence, one can write

$\mathbb{L}^{\star}(\ket{\psi_{0}})=G_{\ket{\psi_{0}}}.$ (19)

We now state our main result.

######  Theorem 1.

Let
$U(t)=\mathcal{O}_{op}\big{[}\prod_{j=1}^{p}\prod_{i=1}^{m}U_{i}(\tau_{j})\big{]}$
be a quantum simulation unitary of duration $t$ as in Eq. (2), divided into
$p$ time-steps, where each $U_{i}(\tau_{j})$ is $k$-local and is time-step
parameterized by $\tau_{j}$, and some operator ordering (such as time-
ordering) is specified. Let $\mathcal{L}_{U_{i}}$ be the string edit map
corresponding to any available parameterization of $U_{i}(\tau_{j})$ and the
set of such maps
$\mathbb{L}\equiv\\{\mathcal{L}_{U_{i}}\>:\>i\in\\{1,\ldots,m\\}\\}$. Let
$B(\mathcal{H}_{2^{n}})\equiv\\{\ket{b}\\}$ be the basis in which computations
(measurements) are being performed classically (quantumly). Then, given an
input state $\ket{\psi_{0}}$, expressed in the basis $B$, if $\ket{b}\notin
G_{\ket{\psi_{0}}}$, where $G_{\ket{\psi_{0}}}$ is constructed according to
Eqs. (17)-(18), then $\bra{b}U(t)\ket{\psi_{0}}=0$.

For a proof, see Appendix B. Note that for $\bra{b_{f}}U(t)\ket{\psi_{0}}$ to
vanish under these conditions, $G_{\ket{\psi_{0}}}$ must satisfy Def. 1 as
demonstrated in Lem. 0.1. In other words, for $G_{\ket{\psi_{0}}}$ to be able
to exclude particular basis states for arbitrary evolution times, it must be a
symmetry-protected subspace. Note that Thm. 1 immediately provides two
corollaries.

Corollary 1 (Post-Selection). For simulation on an idealized, noise-free
quantum computer, if $\ket{b_{f}}\notin G_{\ket{\psi_{0}}}$ then
$||\bra{b_{f}}U(t)\ket{\psi_{0}}||^{2}=0$. Hence, if the state $\ket{b_{f}}$
is measured in the output of a noisy quantum device, it can be assumed that
the state arose as a result of error, and may be discarded.

Corollary 2 (Global Subspace). We assumed a Trotterized form for $U(t)$ in the
statement and proof of Theorem 1. However, we can formally recover the
corresponding global simulation unitary by taking the limit
$p\rightarrow\infty$ where $\tau_{j}=jt/p$ in the time-ordered case and
$\tau_{j}=t/p\>\forall j$ when time-ordering is unnecessary (such as when the
Hamiltonian is time-independent). Nothing in the proof of Theorem 1 relies
upon the finiteness of $p$ or discreteness of the corresponding time
differential $t/p$. Therefore, our result holds for global simulation
unitaries as well. This implies that one can reduce the resource requirements
in classical simulations of $U(t)$ by only evolving basis states
$\ket{b_{f}}\in G_{\ket{\psi_{0}}}$.

###  IV.2 Pseudocode

With an understanding of the recurrence relation in Eq. (17) and how it can
compute symmetry-protected subspaces, we present an algorithm which can
enumerate these subspaces using a breadth-first search. Breadth-first search
to enumerate an entire graph $(V,E)$ has computational complexity
$O(|V|+|E|)$. In our implementation, there are $|G_{\ket{\psi_{0}}}|$ vertices
and we check for $m$ edges at each vertex, giving
$O(m\times|G_{\ket{\psi_{0}}}|)$ edges in the entire graph. Thus, our breadth-
first search to enumerate the symmetry-protected subspace is
$O(|G_{\ket{\psi_{0}}}|+m\times|G_{\ket{\psi_{0}}}|)=O((m+1)\times|G_{\ket{\psi_{0}}}|)$

Algorithm 1 Enumerate symmetry-protected subspace $G_{\ket{\psi_{0}}}$ with a
breadth-first search

String edit maps $\mathbb{L}$ of the simulation operator, initial state
$\ket{\psi_{0}}$

$G\leftarrow\\{\ket{\psi_{0}}\\}$

Let $Q$ be a first-in-first-out queue

$Q\text{.enqueue(}\psi_{0})$

while $Q\neq\emptyset$ do

$\ket{b}\leftarrow Q\text{.dequeue()}$

for
$\ket{b^{\prime}}\in\\{\mathcal{L}_{U_{i}}(\ket{b})\>:\>\mathcal{L}_{U_{i}}\in\mathbb{L}\\}$
do

if $\ket{b^{\prime}}\notin G$ then

$G\leftarrow G\cup\ket{b^{\prime}}$

$Q\text{.enqueue(}\ket{b^{\prime}})$

end if

end for

end while

return $G_{\ket{\psi_{0}}}\leftarrow G$

Algorithm 1 uses the set $G$, which is eventually the symmetry-protected
subspace, to track which states have already been added to $Q$ during the
runtime of the algorithm and prevent them from being checked more than once.
If this set uses the hash of the basis state’s bitstrings, insertion and
search will be average-case $O(1)$. Getting the set of single-operator
transitions $T_{\ket{b}}^{1}$ is $O(m)$, where $m$ is the number of unitary
operators in the system, for a single state $\ket{b}$; while not all $m$
substring edit maps will provide an edge, as many might act as identity on the
state $\ket{b}$, each string edit map must still be checked. This set is
computed for each state discovered, and each state discovered is never added
to the queue $Q$ more than once, which confirms our original complexity
analysis of $O((m+1)\times|G_{\ket{\psi_{0}}}|)$.

This algorithm computes a set equivalent to that described by Eq. (17). See
Appendix B for a proof that this set is a symmetry-protected subspace.

FIG. 3: The number of states in symmetry-protected subspaces for each model,
on a $\log_{2}$ scale. Each blue dot is the size of an individual subspace.
(a) Heisenberg-XXX. (b) $T_{6}$ QCA. (c) $F_{4}$ QCA, which has no known
symmetry operator.

###  IV.3 Usage and limits

As stated in Sec. I, post-selection for quantum simulations is the aim of our
methods. In order to perform post-selection with the algorithm from this
section, the entire symmetry-protected subspace $G_{\ket{\psi_{0}}}$ must be
known, and each measurement $\ket{b_{f}}$ is verified via

$\begin{cases}\ket{b_{f}}\in G_{\ket{\psi_{0}}}&\text{assumed no error}\\\
\ket{b_{f}}\notin G_{\ket{\psi_{0}}}&\text{known error.}\end{cases}$ (20)

Because this requires the computation of the entire subspace, it can still
become computationally intractable. As mentioned in Sec. I, we benchmark this
algorithm against the Heisenberg-XXX, T6 QCA, and F4 QCA quantum simulations.
The computationally limiting factor of this algorithm is the size of the
symmetry-protected subspace, as it has time complexity
$O((m+1)\times|G_{\ket{\psi_{0}}}|)$ and spatial complexity
$O(|G_{\ket{\psi_{0}}}|)$. See Fig. 3 for a depiction of the size of the
subspaces, up to 17 qubits. Figure 3(c) is especially significant because, as
alluded to in Sec. II.2.3, this model previously had no known conservation
laws, but this data shows the partitioning of Hilbert space into symmetry-
protected subspaces. By examining Fig. 3, we can see that the worst case of
each subspace size, $\max(\\{|G|\\})$, is
$\log_{2}(\max(\\{|G|\\}))\approx\log_{2}(|\mathcal{H}_{d}|)-k$ when the model
is comprised of $k$-local operations. The important thing to see is that the
size of symmetry-protected subspaces still scales exponentially in the worst
and average cases; despite them being smaller than the full Hilbert space,
they are only linearly smaller. We address this stop to post-selection with
another algorithm in Sec. V.

##  V Algorithm: Verification of a Shared symmetry-protected Subspace

Because the worst-case size of a symmetry-protected subspace is still
exponential in the number of qubits, Alg. 1 is only applicable at relevant
system sizes for simulations where $G_{\ket{\psi_{0}}}$ is not exponentially
large. Thus, to generalize the usability of symmetry-protected subspaces to
any simulation, we present an alternative algorithm in this section that uses
an efficient, but greedy, heuristic. Instead of computing the entire
exponentially large SPS, this method performs an efficient search for a path
of substring edit maps to connect the initial state $\ket{\psi_{0}}$ and a
measured state $\ket{b_{f}}$ in the graph $D_{\mathbb{L}}$. The heuristic
nature of this algorithm means that it produces only approximate results
(albeit with a degree of approximation that can be continuously improved at
the cost of more time complexity), and since its runtime scales favorably, it
can be used to check the output of quantum simulations well beyond the scale
where other classical methods become intractable.

Naive search algorithms in a graph traditionally use a breadth-first search
from an initial vertex [26], which is what we described in Alg. 1 to enumerate
the symmetry-protected subspace through transitive closure; as stated, this is
too computationally expensive in many cases.

Our greedy algorithm works as follows: Following all edges (application of the
set of substring edit maps $\mathbb{L}$) from a state generates a set of
possible states accessible from that state. Working both forward from
$\ket{\psi_{0}}$ and backward from $\ket{b_{f}}$, we have found a path when
steps from both directions lead to a shared element. To detect whether this
has happened, we use a simple observation about ordered sets: two sets are the
same set if they have the same minimal (or maximal) element. A natural order
for sets of quantum computational basis states (i.e., binary strings) is just
the integer they encode. So, our greedy algorithm works by building two sets
(one starting from the initial state and one from the measured state) and
comparing them; these sets are subsets of the symmetry-protected subspace
corresponding to the initial and final states. If these sets have any common
elements (which we check in constant time by looking at their minimal
elements), then we have found a path connecting the initial and measured
state, and they occupy the same SPS.

In our methods, post-selection’s accept/reject decision is made by checking
for a measurement result in the symmetry-protected subspace of the initial
state. Using the rationale outlined above, we can declare with certainty that
two states share a SPS when their paths collide. On the other hand, if the two
paths terminate in dissimilar minima we assume the initial and measured states
inhabit disjoint subspaces. We say “assume” because the accuracy of our
heuristic differs depending on the simulation; when the locally-minimal choice
at each step from either state does not build a path to the true minimum
element, the two states may build distinct paths while occupying the same SPS.

###  V.1 Searching in the string interaction graph

This algorithm builds a single branch, following local minima in depth-limited
breadth-first searches, from a single starting vertex in the $D_{\mathbb{L}}$
graph, where each vertex $\ket{b}$’s edgeset is given by
$\mathbb{L}(\ket{b})$.

Throughout this section we use the integers encoded by the bitstrings $b$ of
basis states $\ket{b}\in B(\mathcal{H}_{2^{n}})$, which provides a natural
ordering to states. Let the function $\min(A)$ on a set of states $A$ return
the state with the smallest encoded integer in that set:

$\min(A)=\ket{b}\>s.t.\>b\leq b^{\prime},\,\forall\,\ket{b^{\prime}}\in A.$
(21)

To find the minimal element of the symmetry protected subspace, we start with
an element $\ket{b_{0}}=\ket{\psi_{0}}$ or $\ket{b_{f}}$ and build a set
towards the minimum of $G_{\ket{b_{0}}}$, notated $\min(G_{\ket{b_{0}}})$,
with locally optimal decisions. At each step $\ket{b_{\text{curr}}}$ in the
search, we compute $T_{\ket{b_{\text{curr}}}}^{\mu}$ from Eq. (17), which is a
breadth-first search to depth $\mu$, or every state in
$\mathbb{L}^{\mu}(\ket{b_{\text{curr}}})$. Then, the starting point for the
next step, $\ket{b_{\text{next}}}$, is the state with the smallest binary
encoded integer in $T_{\ket{b_{\text{curr}}}}^{\mu}$, notated
$\min(T_{\ket{b_{\text{curr}}}}^{\mu})$. Repeat until the set
$T_{\ket{b_{\text{curr}}}}^{\mu}$ does not offer a state smaller than
$\ket{b_{\text{curr}}}$. Let this process be represented by the recursive
function $\chi$:

$\chi(\ket{b_{0}},\mu)=\begin{cases}\ket{b_{\text{curr}}}&\text{if}\ket{b_{\text{curr}}}=\min(T_{\ket{b_{\text{curr}}}}^{\mu})\\\
\chi(\min(T_{\ket{b_{\text{curr}}}}^{\mu}),\mu)&\text{otherwise}\end{cases}.$
(22)

To reiterate, because each step is given by applications of the substring edit
maps, $\chi(\ket{b_{0}},\mu)$ and $\ket{b_{0}}$ must share a symmetry-
protected subspace, i.e., $\chi(\ket{b_{0}},\mu)\in G_{\ket{b_{0}}}$.
Therefore, when $\chi(\ket{\psi_{0}},\mu)=\chi(\ket{b_{f}},\mu)$ there must be
a sequence of substring edit maps between $\ket{\psi_{0}}$ and $\ket{b_{f}}$;
i.e., a there is a path
$\\{\ket{\psi_{0}}\leftrightarrow\dots\leftrightarrow\ket{b_{f}}\\}\subseteq
D_{\mathbb{L}}$ and $\ket{b_{f}}\in\mathbb{L}^{*}(\ket{\psi_{0}})$.

When both searches conclude in the true minimal element of their corresponding
symmetry-protected subspaces,
$\chi(\ket{\psi_{0}},\mu)=\min(G_{\ket{\psi_{0}}})$ and
$\chi(\ket{b_{f}},\mu)=\min(G_{\ket{b_{f}}})$, we can conclude with certainty
that the states do or do not inhabit the same subspace if
$\chi(\ket{\psi_{0}},\mu)=\chi(\ket{b_{f}},\mu)$ or
$\chi(\ket{\psi_{0}},\mu)\neq\chi(\ket{b_{f}},\mu)$. However, if either search
does not conclude in their targeted minimal state, our conclusions can be
wrong. Suppose $\ket{b_{f}}\in G_{\ket{\psi_{0}}}$, which means that
$\min(G_{\ket{\psi_{0}}})=\min(G_{\ket{b_{f}}})$, but the search result
$\chi(\ket{b_{f}},\mu)$ finds a false minima of $G_{\ket{b_{f}}}$, meaning
$\chi(\ket{b_{f}},\mu)\neq\min(G_{\ket{b_{f}}})$. When post-selecting with our
procedure under these conditions, state $\ket{b_{f}}$ would be wrongly
rejected because the underlying assumption for this heuristic, that
$\chi(\ket{b_{0}},\mu)=\min(G_{\ket{b_{0}}})$, is wrong. In practice, the true
minimal element of the symmetry-protected subspace is unknown; thus, the
assertion that $\chi$ finds the minimum is always an assumption.

Consequently, if the result of the two searches is a collision,
$\chi(\ket{\psi_{0}},\mu)=\chi(\ket{b_{f}},\mu)$, we know the two states must
share a symmetry-protected subspace, even if the searches are at a false
minima. If the two searches do not find a common element, we assume the
initial and measured states occupy separate symmetry-protected subspaces. As
such, measurement outcomes that lie within the intial state’s symmetry-
protected subspace can be rejected. Thus, each measurement $\ket{b_{f}}$ is
verified with

$\begin{cases}\chi(\ket{b_{f}},\mu)=\chi(\ket{\psi_{0}},\mu)&\text{assumed no
error}\\\ \chi(\ket{b_{f}},\mu)\neq\chi(\ket{\psi_{0}},\mu)&\text{assumed
error}\end{cases}.$ (23)

Formally, we can state that the protected subspace formed by $\chi$ around an
initial state $\ket{\psi_{0}}$ is

$P_{G}^{\prime}=\sum_{b}\ket{b}\bra{b}\;s.t.\;\chi(\ket{b},\mu)=\chi(\ket{\psi_{0}},\mu)$
(24)

which is approximately equal to the true symmetry-protected subspace,
$P_{G}\approx P_{G}^{\prime}$.

The confidence of this assumption depends on the quantum system being studied,
and the depth $\mu$ of the local breadth-first searches. Thus, we present
arguments to support it for our three exemplar systems: Heisenberg-XXX,
$T_{6}$ QCA, and $F_{4}$ QCA as outlined in Sec. II.2. For the Heisenberg-XXX
model this assumption is always correct at $\mu=1$: systems whose substring
edit maps $\mathbb{L}$ are isomorphic to nearest-neighbor SWAP substring edit
maps $\mathbb{L}_{\text{SWAP}}$, such as the Heisenberg-XXX model, provably
find $\min(G_{\ket{b_{0}}})$ using $\chi(\ket{b_{0}},\mu=1)$, as shown in
Appendix C.1. We also numerically find that $\chi(\ket{b_{0}},\mu=2)$ is exact
for the $T_{6}$ QCA model. Out of our three studied systems the $F_{4}$ QCA
demonstrates the worst accuracy; $\chi(\ket{b_{0}},\mu=9)$ still returns false
minima, which potentially causes a false rejection. In Sec. V.3.1 we analyze
these properties of the $T_{6}$ and $F_{4}$ QCAs more closely. In Sec. VI.3 we
show that post-selection is still effective when $\chi$ causes a false
rejection to occur.

The computational performance of this method is also model-dependent. Let
$|Depth|$ be the depth of our search, or the number of intermediate states
traversed to (the second line in Eq. (22)) before finding a locally optimal
vertex. This means the sequence of substring edit maps between $\ket{b_{0}}$
and $\chi(\ket{b_{0}},\mu)$ is length $O(\mu\times|Depth|)$, and the
computational complexity of Eq. (22) is $O(m^{\mu}\times|Depth|)$.

Because the reliability and computational complexity of this post-selection
method varies, before applying the algorithm to a given quantum simulation,
benchmarks should be performed on a small version of the system to understand
the $\mu$ required to obtain reliable results.

###  V.2 Pseudocode

Here we describe the algorithm, Alg. 1, to find the minimum binary encoded
state in the symmetry-protected subspace of a given basis state. This
algorithm starts at a state $\ket{b_{0}}$ and follows the local minima of
limited breadth-first searches in an attempt to find the set’s minimal
element, resulting in $\chi(\ket{b_{0}},\mu)$.

Notice that the innermost while loop is borrowed from Alg. 1 to compute each
$T_{\ket{b_{\text{curr}}}}^{\mu}$, where $\ket{b_{\text{curr}}}$ is the
current best minimum state. This time, the queue $Q$ uses the ordered tuple
$(\ket{b^{\prime}},\eta)$ as its elements, where $\ket{b^{\prime}}$ is the
current, un-checked, state in the breadth-first search and $\eta$ is the
number of edges (i.e., depth) that $\ket{b^{\prime}}$ is away from the current
optimal solution $\ket{b_{\text{curr}}}$. Each node checked in this loop has
$O(m)$ new edges, and we check a maximum depth of $\mu$, giving the inner loop
computational complexity $O(m^{\mu})$. The set
$T^{\mu}_{\ket{b_{\text{curr}}}}$ is checked for every node traversed to in
the search, and our search traverses to $|Depth|$ elements, which is how we
arrive at an overall complexity of $O(m^{\mu}\times|Depth|)$ to compute
$\chi(\ket{b_{0}},\mu)$.

Algorithm 1 Greedy pathfinding to the minimum binary encoded state

Path starting state $\ket{b_{0}}$, substring edit maps $\mathbb{L}$, search
limiting integer $\mu$

$\ket{b_{\text{next}}}\leftarrow\ket{b_{0}},\,\ket{b_{\text{curr}}}\leftarrow$
null

while $\ket{b_{\text{curr}}}\neq\ket{b_{\text{next}}}$ do

$T_{\ket{b_{\text{curr}}}}\leftarrow\\{\ket{b_{\text{curr}}},\ket{b_{\text{next}}}\\}$

$\ket{b_{\text{curr}}}\leftarrow\ket{b_{\text{next}}}$

Let $Q$ be a queue

$Q\text{.enqueue}(\ket{b_{\text{curr}}},0)$ $\triangleright$ Each element in
$Q$ is tuple representing (state, depth)

while $Q\neq\emptyset$ do $\triangleright$ Breadth-first search centered
around $\ket{b_{\text{curr}}}$

$\ket{b^{\prime}},\eta\leftarrow Q\text{.dequeue()}$

if $\eta=\mu$ then break $\triangleright$ Reached the depth limit

else

for
$\ket{b^{\prime\prime}}\in\\{\mathcal{L}_{U_{i}}(\ket{b^{\prime}})\>:\>\mathcal{L}_{U_{i}}\in\mathbb{L}\\}$
do

if $\ket{b^{\prime\prime}}\notin T_{\ket{b_{\text{curr}}}}$ then

$Q\text{.enqueue}(\ket{b^{\prime\prime}},\eta+1)$

$T_{\ket{b_{\text{curr}}}}\leftarrow
T_{\ket{b_{\text{curr}}}}\cup\ket{b^{\prime\prime}}$

end if

end for

end if

end while

$\ket{b_{\text{next}}}\leftarrow\min(T_{\ket{b_{\text{curr}}}})$

end while

return $\chi(\ket{b_{0}},\mu)=\ket{b_{\text{curr}}}$

To post-select with Alg. 1, run it once for $\ket{\psi_{0}}$ to get
$\chi(\ket{\psi_{0}},\mu)$, again for each measurement $\ket{b_{f}}$ to get
$\chi(\ket{b_{f}},\mu)$, and check the results in Eq. (23).

As an aside, if all states computed in the search,
$T^{\mu}_{\ket{b_{\text{curr}}}}$, are cached with the result of their search,
subsequent executions of the algorithm can be preempted with their previously
calculated result if a state is found inside the cache.

###  V.3 Benchmarking data for Algorithm 1

This subsection presents benchmark data for Alg. 1 that demonstrates its
practical applicability for our example systems. These benchmarks include
heuristic reliability and search depth. We supplement this data with proofs or
explicit analytical forms when possible. When this is not possible, we instead
rely on extrapolation from the data to inform the asymptotics of these
quantities. We find that for the Heisenberg-XXX model, the heuristic is
provably exact (Appendix C.1) and the search depth has an explicit equation
(Appendix C.2), while the $T_{6}$ and $F_{4}$ QCA models rely on data-driven
intuition for these quantities.

FIG. 4: The proportion of states $\ket{b}$ where
$\chi(\ket{b_{0}},\mu)\neq\min(G_{\ket{b}})$ via Alg. 1, called “failed
searches”, at different system sizes, for each $\mu$ used. (a) $T_{6}$ QCA. At
$\mu=2$, Alg. 1 becomes exact for this system. (b) $F_{4}$ QCA. High $\mu$ can
still fail, but does asymptotically better for $\mu\geq 5$.

####  V.3.1 Search reliability

If searching for the minimal element from $\ket{b_{0}}=\ket{\psi_{0}}$ or
$\ket{b_{f}}$ fails, because $\chi(\ket{b_{0}},\mu)\neq\min(G_{\ket{b_{0}}})$,
then it is possible to falsely assume that two states occupy disjoint
subspaces. As such, we examine the reliability of the searching function
$\chi$ on each exemplar system.

In the Heisenberg-XXX model, the base case for $\chi$, as shown in Eq. (22),
only activates at a state which is always the minimum binary-encoded state in
the symmetry-protected subspace. This is because the set of substring edit
maps for the Heisenberg-XXX system is isomorphic to the set of substring edit
maps for the nearest-neighbor SWAP network:
$\mathbb{L}_{\text{HeisXXX}}=\mathbb{L}_{\text{SWAP}}$. Appendix C.1 proves
that $\chi$ at $\mu=1$ is exact for $\mathbb{L}_{\text{SWAP}}$.

Because the $T_{6}$ QCA and $F_{4}$ QCA do not have SWAP-isomorphic substring
edit maps, we do not have any available proofs for their exactness, and
instead rely on simulating $\chi$ for these models. The results of this are in
Fig. 4, where for the $T_{6}$ and $F_{4}$ QCA models we plot the proportion of
states $\ket{b}\in B(\mathcal{H}_{d})$ such that
$\chi(\ket{b},\mu)\neq\min(G_{\ket{b}})$ for each $n\in[k,17]$.

First, notice Fig. 4(a), which shows search failures for the $T_{6}$ QCA
model; while $\mu=1$ causes most searches to fail as $n$ increases, $\mu=2$
causes every search to succeed for every $n$ in the domain. Thus, it we can
assume that the heuristic for Alg. 1 is accurate for this model and that
$\chi(\ket{b},\mu=2)=\min(G_{\ket{b}}),$ $\forall\ket{b}\in B$.

Next, notice Fig. 4(b), which shows search failures for the $F_{4}$ QCA model;
we use $\mu=1,3,5,7,9$, and find that $\chi$ can fail even with high $\mu$. We
see that $\mu=1,3$ sees majority failures over the domain, that $\mu=5,7,9$
has bounded failure $\approx 0.05$, and that no $\mu$ completely removes
failures.

This metric only determines the reliability of $\chi$, not the reliability of
post-selecting with $\chi$. This is because, as discussed previously in Sec.
V.1, two searches can fail and obtain the same false minima,
$\chi(\ket{\psi_{0}},\mu)=\chi(\ket{b_{f}},\mu)\neq\min(G_{\ket{\psi_{0}}})=\min(G_{\ket{b_{f}}})$,
i.e., two failed searches can still successfully identify that the states
share a conserved quantity. Therefore, the search reliability shown in Fig. 4
should be seen as an upper-bound on the error caused by the heuristic.

Sections VI.2 and VI.3 explore this reliability further in the context of
post-selection on noisy quantum data. In the context of post-selection, we
find that even small $\mu=3$ is acceptable for the $F_{4}$ QCA, as the most
likely measurement outcomes also share false minima with the initial state. We
discuss this observation more with the data that supports it in Sec. VI.3.

FIG. 5: The longest path taken at each $n$ for the models (a) Heisenberg-XXX
(b) $T_{6}$ QCA (c) $F_{4}$ QCA.

####  V.3.2 Search depth

The time performance of Alg. 1 is heavily dependent on the length of the path
taken. The asymptotic worst-case of the search depth is system-dependent,
being identified as $n^{2}/4$ for the Heisenberg-XXX model in Appendix C.2,
and closer to $O(n\ln(n))$ for the $T_{6}$ QCA and $O(n\log_{10}(n))$ for the
$F_{4}$ QCA by observing the largest search depth for the algorithm at small
$n$, shown in Fig. 5. In this figure we find the largest path length at each
system size for each relevant $\mu$. For the $T_{6}$ model, shown in Fig.
5(a), search depth becomes longer as the parameter $\mu$ increases, indicating
that false minima are avoided. For the $F_{4}$ model, shown in Fig. 5(b), the
search depth becomes shorter as $\mu$ increases; this is because higher $\mu$
means more edges are traversed at each step in the search, and thus less depth
is needed in the search overall.

These path lengths, combined with the $\mu$ found in Sec. V.3.1, results in
the computational complexity of Alg. 1 on a single state being: Heisenberg-XXX
is $O(mn^{2})$, $T_{6}$ QCA is $O(m^{2}n\ln(n))$, and $F_{4}$ QCA is
$O(m^{\mu}n\log_{10}(n))$. The $\mu$ exponent is left in the $F_{4}$ QCA
because, as shown in Sec. V.3.1, we find no value of $\mu$ which guarantees
that a search finds the minimal element, and thus expect $\mu$ to vary
depending on the symmetry-protected subspace occupied for the simulation. As
we can see in Fig. 4(b), $\mu\geq 5$ is mostly sufficient.

FIG. 6: The fidelities $\mathcal{F}(P_{\text{sim}}(p),P_{\text{ideal}}(p))$,
$P_{\text{sim}}=P_{\text{raw}},P_{\text{ps},\mu=1},P_{\text{ps},\mu=2},P_{\text{ISPS}}$
of each model (rows) at errors $\epsilon/3=0.005,0.01,0.02,0.05$ (columns).
There are $10,000$ measurements at each Trotter step. The $\mu$ used is the
smallest required to reach perfect shared subspace verification under Alg. 1.
This figure shows that with post-selection, model fidelity remains above
random noise at least as long as $p=2*n$.

##  VI Post-Selection in Emulated Noisy Simulations

To demonstrate the power of our algorithm, we perform automated post-selection
on simulated quantum computations with noise. We perform discretized time
evolutions with our exemplar systems, the Heisenberg-XXX, $T_{6}$ QCA, and
$F_{4}$ QCA models, and show that our post-selection methods restructure noisy
measurement probability distributions to be closer to the one found in an
ideal simulation. We first run the simulation without errors via
$U(t)\ket{\psi_{0}}$, constructing the measurement distribution at each time
step. Next, we repeated the simulations with depolarizing noise injected after
each layer of gates with independent probability
$\epsilon/3=0.005,0.01,0.02,0.05$ of $X_{i},Y_{i},Z_{i}$ on each qubit $i$.
See Appendix D.1 for more noise model details. We use Kullback-Liebler
divergence [27] to quantify the distance between the ideal and noisy
measurement distributions constructed at each measurement layer, with and
without post-selection using Alg. 1 and the symmetry-check in Eq. (23).

We observe a significant increase in accuracy of the data when using our post-
selection methods without adding circuit runs, even in simulation subspaces
which have imperfect pathfinding.

###  VI.1 Methods

We simulate discrete time evolution, and for each discrete time step a
sequence of measurements is used to construct a probability distribution. Let
$P(p)$ be the measurement distribution of the wavefunction after $p$ Trotter
steps, such that $P(p)\approx||\ket{\psi(p)}||^{2}$ up to shot noise. We
compare an ideal simulation measurement distribution $P_{\text{ideal}}(p)$, a
noisy simulation without post-selection $P_{\text{raw}}(p)$, and a noisy
simulation with post-selection using Alg. 1 at a given $\mu$,
$P_{\text{ps},\mu}(p)$. Let $P(p)$ be a normalized probability distribution
constructed by a sequence of $M$ measurements after $p$ Trotter steps, and
$P(p,b)\approx||\bra{b}\ket{\psi(p)}||^{2}$, which is the probability
amplitude of state $\ket{b}$ at that layer $p$, up to shot noise. We use
Kullback-Liebler divergence [27],

$\mathcal{D}(P(p),Q(p))\equiv\sum_{b\in
B(\mathcal{H}_{d})}P(p,b)\ln(\frac{P(p,b)}{Q(p,b)}),$ (25)

which measures the distance between two probability distributions, and use the
following equation from [28] as a simulation’s fidelity,

$\mathcal{F}(P_{\text{sim}}(p),P_{\text{ideal}}(p))\equiv
1-\frac{\mathcal{D}(P_{\text{sim}}(p),P_{\text{ideal}}(p))}{\mathcal{D}(P_{\text{IRN}},P_{\text{ideal}}(p))},$
(26)

where $P_{\text{IRN}}$ is the incoherent random noise probability distribution

$P_{\text{IRN}}(b)\equiv\frac{1}{d}\>\>\forall\>\>\ket{b}\in
B(\mathcal{H}_{d})$ (27)

which we expect to measure once noise has proliferated in the computation.
Equation (26) will return 1 if $P_{\text{sim}}=P_{\text{ideal}}$, and 0 if
$P_{\text{sim}}=P_{\text{IRN}}$. We compute Eq. (26) at each measurement layer
in the computation, with $P_{\text{sim}}$ being either the raw data
$P_{\text{raw}}$ or the post-selected data $P_{\text{ps}}$.

We post-select by removing measurement results that fail under Alg. 1 and Eq.
(23) and re-normalizing the probability distribution. If $P_{\text{raw}}(p)$
is the measurement distribution with noise at step $p$ of $M$ measurements,
$P_{\text{ps}}(p)$ is the post-selected measurement distribution, given by
checking each state in $P_{\text{raw}}(p)$ in Alg. 1, at the same Trotter-step
$p$ with measurements $M_{\text{kept}}\leq M$.

FIG. 7: The maximum, minimum, mean, and standard deviation of the set of
fidelities, $\\{\mathcal{F}\\}$, generated by simulating the $F_{4}$ QCA for
each initial condition in its rocky subspaces at $\mu=9$ with post-selected
distributions
$P_{\text{sim}}=P_{\text{ps},\mu=5},P_{\text{ps},\mu=7},P_{\text{ps},\mu=9},P_{\text{SPS}}$.
The error rate is (a) $\epsilon/3=0.02$ and (b) $\epsilon/3=0.0$, or an ideal
simulation.

###  VI.2 Post-selection with perfect symmetry-protected subspaces

We run a 15 qubit simulation of each of the three exemplar models, using
initial conditions seen in the literature, the anti-ferromagnetic state
$\ket{\psi_{0}}=\ket{1010\ldots 101}$ for the Heisenberg-XXX model, a single
bit-flip on the middle qubit $\ket{\psi_{0}}=\ket{0\ldots 1\ldots 0}$ for the
$T_{6}$ QCA [21, 22], and two bit-flips neighboring the center qubit
$\ket{\psi_{0}}=\ket{0\ldots 101\ldots 0}$ for the $F_{4}$ QCA [22]. No false
minima are seen past $\mu=2$ in any of the simulations for these initial
conditions, meaning post-selection functions effectively. The data uses
$M=10,000$ measurements at each measurement layer and does not add any
measurements to replace error-victim circuit runs.

Results can be seen in Fig. 6. We observe that post-selected data approaches
incoherent noise in the SPS of that initial state $P_{\text{ISPS}}$, defined
as

$P_{\text{ISPS}}(b)\equiv\begin{cases}\frac{1}{|G|}&\>\ket{b}\in G\\\
0&\>\ket{b}\notin G\\\ \end{cases}$ (28)

Because our simulation’s raw data converges to random noise, which is a
uniform probability distribution over the Hilbert space, and post-selection
only removes measurement outcomes, we should expect that the result of post-
selection is still a uniform probability distribution, just over the symmetry-
protected subspace instead. In other words, we should observe that

$\lim_{p\rightarrow\infty}\mathcal{F}(P_{\text{ps}}(p),P_{\text{ideal}}(p))\approx\mathcal{F}(P_{\text{ISPS}}(p),P_{\text{ideal}}(p))$
(29)

given enough measurements.

###  VI.3 Post-selection with incorrect symmetry-protected subspaces

We show that even in subspaces of the $F_{4}$ QCA model where Alg. 1 has
failures at high $\mu$, which is equivalent to false-positive mislabeling some
basis vectors as outside the symmetry-protected subspace, our methods mostly
mitigate more errors than they introduce and thus still show promise for NISQ
quantum simulations.

We found 18 symmetry-protected subspaces of the $F_{4}$ QCA that had false
minima at $n=15$ with $\mu\geq 9$ (henceforth called rocky subspaces) and ran
simulations to $p=29$ with $M=5000$ measurements at each cycle, initialized in
one state from each rocky subspace. See Appendix D.4.2 for a list of these
initial states. This way, Alg. 1 is highly likely to meet false minima and
reject measurements that were in the symmetry-protected subspace of the
initial state. We calculate
$\mathcal{F}(P_{\text{ps},\mu}(p),P_{\text{ideal}}(p))$ for each simulation,
and run our post-selection method with $\mu\in=5,7,9$. See Fig. 7(a) for the
set of fidelities $\\{\mathcal{F}\\}$ at $\epsilon/3=0.02$, compared to post-
selection with the “Full SPS”, which is post-selection using the full
symmetry-protected subspace, instead of the heuristic; in the average case our
post-selection keeps the simulation data above incoherent noise, even with
$\mu=5,7$. The worst-case $\mathcal{F}$ (in red) drops below the statistical
difference of incoherent noise, which indicates that our methods would not
work on this particular simulation. However, the mean (orange) and standard
deviation (orange bars) remain reliably higher than incoherent noise, with the
maximum $\mathcal{F}$ keeping close to the ideal simulation.

The other important observation is that our post-selection, which is using an
inaccurate symmetry-protected subspace, obtains a higher mean and maximum
fidelity than post-selection with a well-defined symmetry-protected subspace.
We theorize the reason for this is as follows: if two states share a false
minima with $\chi$, they are likely close to each other in the string
interaction graph $D_{\mathbb{L}}$, and thus have a higher probability of
transitions between each other in the simulation. On the other hand, if a
measured state does not share a false minima with the initial state, it is
less likely to be observed in a simulation, but a noisy simulation will
artificially amplify the probability of that state. The result is that post-
selection with the true SPS will accept these artificially amplified results,
where our false SPS will coincidentally reject states that can possibly have
this happen.

We also examine our post-selection in a simulation absent any error; this is
to see how much a false rejection can degrade a perfect computation. See Fig.
7(b) for post-selection with an error-free simulation. It can be seen that
most fidelities remain above $\approx 0.95$, meaning most of the time this
error-prone post-selection method will have minimal impact even on an ideal
simulation. The worst-case degraded fidelity is $\approx 0.85$.

##  VII Discussion and conclusion

We introduce a graph theory interpretation of quantum time evolution, which
provides a theoretical framework through which symmetry-protected subspaces
can be constructed via transitive closure. We identify that these invariant
subspaces are an operator-free method for characterizing symmetry in the
system, by indicating a conserved quantity of an initial state as it is
manipulated by the dynamical system without needing to explicitly identify it.

We observe that a symmetry-protected subspace can be used to provide a smaller
computational space, post-select noisy quantum simulation data, or be analyzed
to deduce a symmetry operator. We identify post-selection as a pertinent
application and introduce two main classical algorithms that elucidate the
features of a quantum system’s symmetry-protected subspaces. These algorithms
employ a basis string edit map, which is an efficient construction to provide
the local dynamics of an operator by focusing on the presence or absence of
basis vectors through its action.

The first algorithm uses transitive closure, calculated with breadth-first
search of these basis string edit maps, to enumerate every state within the
symmetry-protected subspace of the initial state. Because post-selection with
this subspace requires totalling the entire SPS, which can still be
exponentially difficult, we provide a second, polynomial-scaling algorithm.
The second algorithm attempts to find the smallest basis state (by binary-
encoded integer) in an SPS by following a locally-optimal heuristic, and
establishing that a path of basis string edit maps exists between an initial
and final state if their respective searches collide at any elements.

We conclude by demonstrating post-selection using our second algorithm, which
shows that even when the raw simulation data degrades to incoherent noise, our
methods effectively recover a probability distribution closer to an ideal
computation. Our methods are compatible with any other error-mitigation
technique compatible with post-selection, such as zero noise extrapolation,
and thus present a further addition to a growing array ot techniques to
improve noisy quantum computation [29].

We identify a few obvious extensions to this project. First, if Alg. 1
concludes with certain classical resources, this implies that the wavefunction
can be stored with constant overhead on those same resources; this could
shrink the computational memory requirements from the quantum regime to the
classical regime. Second, we speculate that more reliable algorithms than Alg.
1 may fulfill the same function; a string-matching algorithm similar to
Needleman-Wunsch [30] with the basis string edit maps, instead of
insert/delete/shift edits, may fulfill this function.

This work also presents interesting results on post-selection that are worth
further exploration. The first is the convergence to incoherent noise _within_
a symmetry-protected subspace; this shows that, as expected, post-selection
will not converge to an ideal computation, regardless of the additional
circuit runs to replace error-victim runs detected, as error will still have
proliferated within the the SPS. However, such incoherent distributions still
have significant complexity and nontrivial structure (see e.g., [21]), making
them particularly interesting in cases where the underlying symmetry is not
analytically known. The second is the observation that mis-identifying
measurements within the SPS paradoxically results in a more accurate
computation than what results from full knowledge of the SPS. The exact cause
of this is up to speculation, and is likely a problem-dependent effect, but is
interesting in its own right.

Code and data are available upon reasonable request to the corresponding
authors.

##  VIII Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers,
Superconducting Quantum Materials and Systems Center (SQMS) under contract
number DE-AC02-07CH11359, and by National Science Foundation grant
PHY-1653820. This work was also authored in part by the National Renewable
Energy Laboratory (NREL), operated by Alliance for Sustainable Energy, LLC,
for the U.S. Department of Energy (DOE) under Contract No. DE-AC36- 08GO28308.
This work was supported by the Laboratory Directed Research and Development
(LDRD) Program at NREL. The views expressed in the article do not necessarily
represent the views of the DOE or the U.S. Government. The U.S. Government
retains and the publisher, by accepting the article for publication,
acknowledges that the U.S. Government retains a nonexclusive, paid-up,
irrevocable, worldwide license to publish or reproduce the published form of
this work, or allow others to do so, for U.S. Government purposes.

##  APPENDIX A PROOF OF ACCEPTABLE MAXIMAL COVERAGE IN SYMMETRY-PROTECTED
SUBSPACES

Take the operator $P_{G}=\sum_{b\in G_{\ket{\psi_{0}}}}\ket{b}\bra{b}$ to be
the projection operator for the symmetry-protected subspace of
$\ket{\psi_{0}}$:

$[P_{G_{\ket{\psi_{0}}}},U]=0$ (30)

and suppose it has been constructed such that for each $\ket{b}\in
G_{\ket{\psi_{0}}}$, $\exists\,t\;s.t.\;\bra{b}U^{t}\ket{\psi_{0}}\neq 0$. Now
suppose that the basis vector $\ket{k}$ is never seen by time evolution of
$U$, such that

$\bra{k}U(t)\ket{\psi_{0}}=0\;\forall t.$ (31)

If we include $\ket{k}$ and the symmetry-protected subspace of $\ket{k}$ in
our protected subspace, such that $G^{\prime}=G_{\ket{\psi_{0}}}\cup
G_{\ket{k}}$, then $P_{G^{\prime}}=P_{G_{\ket{\psi_{0}}}}+\sum_{j\in
G_{\ket{k}}}\ket{j}\bra{j}$, and we get

$\displaystyle[P_{G^{\prime}},U]$
$\displaystyle=[P_{G_{\ket{\psi_{0}}}}+P_{G_{\ket{k}}},U]$ (32)
$\displaystyle=[P_{G_{\ket{\psi_{0}}}},U]+[P_{G_{\ket{k}}},U]$ (33)
$\displaystyle=0$ (34)

by the Def. 1 of a symmetry-protected subspace for both $G_{\ket{\psi_{0}}}$
and $G_{\ket{k}}$. This shows that if a state which is never seen by the time
evolution is included in the final symmetry-protected subspace, the projection
operator still commutes and thus the subspace still follows Def. 1.

##  APPENDIX B PROOF OF THEOREM LABEL:THM:SPSTHEOREM

We take $\ket{\psi_{0}}\in B$ to be a product state. Due to linearity,
generalizing the proof to superposition states is trivial, given that one
interprets $G_{\ket{\psi_{0}}}$ as being the union of all symmetry-protected
subspaces of each of which corresponds to one (or more) basis states in the
superposition. For notational simplicity, and without loss of generality, we
also assume an implicit operator ordering $\mathcal{O}_{op}$. Consider the
transition amplitude

$\begin{split}\bra{b_{f}}U(t)\ket{\psi_{0}}&=\bra{b_{f}}\prod_{j=1}^{p}\prod_{i=1}^{m}U_{i}(\tau_{j})\ket{\psi_{0}}\\\
&=\bra{b_{f}}\prod_{j=1}^{p}\prod_{i=1}^{m}\sum_{\ket{b_{j,i}}\in
B}\ket{b_{j,i}}\bra{b_{j,i}}U_{i}(\tau_{j})\ket{\psi_{0}},\end{split}$ (35)

where in the second line of Eq. (35) we have inserted the identity operator
$\mathds{1}=\sum_{\ket{b_{j,i}}\in B}\ket{b_{j,i}}\bra{b_{j,i}}$ after each
$k$-local unitary operator. Without loss of generality, we assume an operator
ordering to expand and rearrange the right-hand side of Eq. (35),

$\begin{split}&\bra{b_{f}}U(t)\ket{\psi_{0}}=\\\ &=\sum_{\ket{b_{p,m}}\in
B}\bra{b_{f}}\ket{b_{p,m}}\sum_{\ket{b_{p,m-1}}\in
B}\bra{b_{p,m}}U_{m}(\tau_{p})\ket{b_{p,m-1}}\ldots\\\
&\ldots\sum_{\ket{b_{1,1}}\in
B}\bra{b_{1,2}}U_{2}(\tau_{1})\ket{b_{1,1}}\bra{b_{1,1}}U_{1}(\tau_{1})\ket{\psi_{0}}.\end{split}$
(36)

Next, we note that by the definition of the string edit map, Def. 4, and the
presence of the local transition amplitude
$\bra{b_{1,1}}U_{1}(\tau_{1})\ket{\psi_{0}}$ in Eq. (36), we can re-write

$\begin{split}&\sum_{\ket{b_{1,1}}\in
B}\bra{b_{1,2}}U_{2}(\tau_{1})\ket{b_{1,1}}\bra{b_{1,1}}U_{1}(\tau_{1})\ket{\psi_{0}}=\\\
=&\sum_{\ket{b_{1,1}}\in\mathcal{L}_{U_{1}}(\ket{\psi_{0}})}\bra{b_{1,2}}U_{2}(\tau_{1})\ket{b_{1,1}}\bra{b_{1,1}}U_{1}(\tau_{1})\ket{\psi_{0}}\\\
=&\sum_{\ket{b_{1,1}}\in
T^{1}_{\ket{\psi_{0}}}}\bra{b_{1,2}}U_{2}(\tau_{1})\ket{b_{1,1}}\bra{b_{1,1}}U_{1}(\tau_{1})\ket{\psi_{0}},\end{split}$
(37)

where in the third line of Eq. (37) we have used the fact that
$\mathcal{L}_{U_{1}}(\ket{\psi_{0}})\subseteq T^{1}_{\ket{\psi_{0}}}\subseteq
B$, which is to say, we can expand the summation to run over
$T^{1}_{\ket{\psi_{0}}}$ by summing over terms where the local transition
amplitude is zero. By the same logic, we can re-write the second summation as

$\begin{split}\sum_{\ket{b_{1,2}}\in
B}(\ldots)&=\sum_{\ket{b_{1,2}}\in\mathcal{L}_{U_{2}}(\ket{b_{1,1}})\>:\>\ket{b_{1,1}}\in\mathcal{L}_{U_{1}}(\ket{\psi_{0}})}(\ldots)\\\
&=\sum_{\ket{b_{1,2}}\in T^{2}_{\ket{\psi_{0}}}}(\ldots),\end{split}$ (38)

where again, we have simplified notation by noting that
$\\{\mathcal{L}_{U_{2}}(\ket{b_{1,1}})\>:\>\ket{b_{1,1}}\in\mathcal{L}_{U_{1}}(\ket{\psi_{0}})\\}\subseteq
T^{2}_{\ket{\psi_{0}}}$ and summing over all additional states in
$T^{2}_{\ket{\psi_{0}}}$ to which the transition amplitude is zero. Generally,
one can re-write summation $(j-1)m+i$ as

$\sum_{\ket{b_{j,i}}\in B}(\ldots)=\sum_{\ket{b_{j,i}}\in
T^{(j-1)m+i}_{\ket{\psi_{0}}}}(\ldots).$ (39)

We therefore write the full transition amplitude as

$\begin{split}&\bra{b_{f}}U(t)\ket{\psi_{0}}=\\\ &\sum_{\ket{b_{p,m}}\in
T_{\ket{\psi_{0}}}^{pm}}\bra{b_{f}}\ket{b_{p,m}}\sum_{\ket{b_{p,m-1}}\in
T_{\ket{\psi_{0}}}^{pm-1}}\bra{b_{p,m}}U_{m}(\tau_{p})\ket{b_{p,m-1}}\\\
&\ldots\sum_{\ket{b_{1,1}}\in
T^{1}_{\ket{\psi_{0}}}}\bra{b_{1,2}}U_{2}(\tau_{1})\ket{b_{1,1}}\bra{b_{1,1}}U_{1}(\tau_{1})\ket{\psi_{0}}.\end{split}$
(40)

There are now two cases to consider. Let $T_{\ket{\psi_{0}}}^{\ell}$ be the
iteration at which the stop condition in Eq. (17) is activated. Either
$T_{\ket{\psi_{0}}}^{pm}=T_{\ket{\psi_{0}}}^{\ell}\equiv G_{\ket{\psi_{0}}}$,
or else $T_{\ket{\psi_{0}}}^{pm}\subset T_{\ket{\psi_{0}}}^{\ell}$. In either
case, $\ket{b_{f}}\notin G_{\ket{\psi_{0}}}\implies\ket{b_{f}}\notin
T_{\ket{\psi_{0}}}^{pm}$, which means that
$\bra{b_{f}}\ket{b_{p,m}}=0\>\>\forall\>\>\ket{b_{p,m}}\in
T_{\ket{\psi_{0}}}^{pm}$ and the full transition amplitude must vanish. ∎

##  APPENDIX C ANALYSIS OF ALGORITHM LABEL:ALG:GREEDYPATH

###  C.1 Proof of exactness for $\mathbb{L}_{\text{SWAP}}$ isomorphic systems

We will show that for a system $U$ which has
$\mathbb{L}_{U}=\mathbb{L}_{\text{SWAP}}$, iterative path creation with Eq.
(22) (which is what Alg. 1 computes) has only one stop condition when using
$\mu=1$, which is the binary encoded minima of the entire symmetry-protected
subspace.

Assume our simulation of $n$ qubits has substring edit maps
$\mathbb{L}_{\text{SWAP}}\equiv\\{\mathcal{L}_{\text{SWAP}_{0,1}},\mathcal{L}_{\text{SWAP}_{1,2}},\ldots,\mathcal{L}_{\text{SWAP}_{n-2,n-1}}\\}$.
We take $\ket{b_{0}}$ to be the first state in the path, and it is a bitstring
$\ket{b_{0}}=\ket{\\{0,1\\}^{n}}$ such that
$\sum_{i=0}^{n-1}a_{i}^{\dagger}a_{i}\ket{b_{0}}=s\ket{b_{0}}$. The smallest
binary encoding belonging to a state in the symmetry-protected subspace of
$\ket{b_{0}}$, which is the SPS encoding particle conservation symmetry, will
belong to the bitstring $\min(G_{\ket{b_{0}}})=\ket{1^{s}0^{n-s}}$. The bits
are ordered left-to-right as least-to-most significant.

The string edit map $\mathcal{L}_{\text{SWAP}_{i,i+1}}$ acting on its local
basis set $B_{i}=\\{\ket{00},\ket{01},\ket{10},\ket{11}\\}$ has the following
mappings:

$\begin{split}\mathcal{L}_{\text{SWAP}_{i,i+1}}(\ket{00})&=\\{\ket{00}\\}\\\
\mathcal{L}_{\text{SWAP}_{i,i+1}}(\ket{01})&=\\{\ket{10},\ket{01}\\}\\\
\mathcal{L}_{\text{SWAP}_{i,i+1}}(\ket{10})&=\\{\ket{01},\ket{10}\\}\\\
\mathcal{L}_{\text{SWAP}_{i,i+1}}(\ket{11})&=\\{\ket{11}\\}\\\ \end{split}.$
(41)

When applied to a $n$-qubit state represented as an $n$-character bitstring
$b$, the mapping operates as e.g.

$\begin{split}\mathcal{L}_{\text{SWAP}_{i,i+1}}(&\ket{b[0,i-1];01;b[i+2,n-1]})\\\
&=\ket{b[0,i-1];10;[i+2,n-1]},\end{split}$ (42)

where the “;” symbol means string concatenation and “$b[i,j]$” is the
substring of $b$ from characters at index $i$ through $j$ (inclusive).

Given an incomplete search with $\chi(\ket{b_{0}},\mu=1)$, where
$\ket{b_{j}}=\min(T^{1}_{\ket{b_{j-1}}})$ is the most recent progression in
the search, the next state $\ket{b_{j+1}}$ in the search is the smallest
bitstring in $\mathbb{L}_{\text{SWAP}}(\ket{b_{j}})$, given by
$\ket{b_{j+1}}=\min(T^{1}_{\ket{b_{j}}})$; this will be
$\mathcal{L}_{\text{SWAP}_{i,i+1}}(\ket{b_{j}})$, where $\ket{b_{j}[i+1]}$ is
the most significant bit equal to $\ket{1}$ such that
$\ket{b_{j}[i]}=\ket{0}$.

In other words, our heuristic dictates that we will only choose
$\mathcal{L}_{\text{SWAP}_{i,i+1}}(\ket{b_{j}})$ when
$\ket{b_{j}[i,i+1]}=\ket{01}$. Therefore, Alg. 1 will continue until there is
no substring “01” in $b_{j}$. For a bitstring, the only configuration of bits
that meets this condition is $b=1^{s}0^{n-s}$, which is the state
$\min(G_{\ket{b_{0}}})$. Therefore, if Alg. 1 halts, it must have found the
minimum binary encoded state for that symmetry-protected subspace. ∎

###  C.2 Max depth of search in $\mathbb{L}_{\text{SWAP}}$

Here we calculate the worst-case search depth for Alg. 1 when using string
edit maps $\mathbb{L}_{\text{SWAP}}$. Each time a new state $\ket{b_{j+1}}$ is
traversed to with $\ket{b_{j+1}}=\min(T^{1}_{\ket{b_{j}}})$ from Eq. (22), a
“1” bit is swapped with a left-neighboring “0” bit, assuming bits are least-
to-most significant left-to-right. We count each of these swaps as 1 step in
the search, and this subsection will show an analytical form for the most
steps ever required for an arbitrary bitstring $\ket{b_{0}}$ of length $n$ to
reach the stop condition $\ket{1^{s}0^{n-s}}$, as defined at the beginning of
Appendix C.1. That section also shows that for the SWAP-isomorphic model we
always get $\chi(\ket{b_{0}},1)=\ket{1^{s}0^{n-s}}$. Note that, because
$\mu=1$ is used, the path length of edges in $D_{\mathbb{L}}$ between
$\ket{b_{0}}$ and $\ket{1^{s}0^{n-s}}$ is equal to the search depth, because
$|Path|=\mu\times|Depth|=|Depth|$.

Suppose we have a symmetry-protected subspace of $s$ many “1” bits and $n-s$
many “0” bits. The target state $\min(G_{\ket{b_{0}}})$ for a search with
$\chi(\ket{b_{0}},1)$ is $\min(G_{\ket{b_{0}}})=\ket{1^{s}0^{n-s}}$. The
longest path will start with $\ket{b_{0}}=\ket{0^{n-s}1^{s}}$, since every
other state is fewer SWAPs away from $\min(G_{\ket{b_{0}}})$.

Each $\ket{1}$ state in $\ket{b_{0}}$ requires $n-s$ swaps to get to its
position in $\ket{1^{s}0^{n-s}}$; therefore, in a system of $s$ “$\ket{1}$”
states, the search depth obeys the summation

$|Depth|=\sum_{i=1}^{s}n-s,$ (43)

which simplifies to

$|Depth|=s(n-s).$ (44)

This quantity is maximized at $s=n/2$, which results in
$\max(|Depth|)=\frac{n^{2}}{4}$.

##  APPENDIX D SIMULATION INFORMATION

Our simulations all follow the same template: we take an initial state
$\ket{\psi_{0}}$, which is a $Z$ basis vector prepared by single-qubit
Pauli-$X$ gates. In our ideal simulations, we follow this with $p\in[0,29]$
gate layers of gates to encode the dynamical system, followed by a parallel-
readout of all qubits in the $Z$ basis with $M$ measurements. This is shown in
the circuit diagram below, where $\ket{\psi_{0}}$ is a computational basis
state prepared by a set of $X$ gates, each $U(\tau_{j})$ is a single layer of
the Trotterized circuit, at layer $j$, with parameters $\tau_{j}$. Once the
system has been simulated up to $t$-layers, we measure $M$ times, and repeat
for each $p\in[0,29]$:

.

For a simulation with noise, we also added an “error layer” $E$ after each
$U(\tau_{j})$, such that $U\rightarrow UE$. It is detailed below in Appendix
D.1.

###  D.1 Simulation noise model

We use symmetric single-qubit depolarizing noise with probability $\epsilon/3$
after each Trotter layer in the circuit. The circuit diagram for this process
at an arbitrary Trotter step $j$ is

,

where a gate $X(\epsilon/3)$ is an $X$ gate applied with probability
$\epsilon/3$ and identity $\mathds{1}$ applied with probability
$1-\epsilon/3$. In practice, we simulate this with density matrices through
the following steps: if we have the wavefunction’s density matrix $\rho$ from
the last application of the Trotterized circuit layer $U(\tau_{j})$, we apply
the symmetric depolarizing noise with the equation

$\rho\rightarrow(1-\epsilon)\rho+\frac{\epsilon}{3}X_{i}\rho
X_{i}+\frac{\epsilon}{3}Y_{i}\rho Y_{i}+\frac{\epsilon}{3}Z_{i}\rho Z_{i}$
(45)

on each qubit $i$. As stated in Sec. VI, we use
$\epsilon/3=0.005,0.01,0.02,0.05$.

###  D.2 Heisenberg-XXX

Here we give the simulation details for the Heisenberg-XXX model, with
Hamiltonian given in Sec. II.2.1 and simulation results in Sec. VI.2.

####  D.2.1 Circuit implementation

The $k$-local unitary used ($k=2$) is
$U_{i,i+1}=\text{iSWAP}_{i,i+1}(\theta)Z_{i}Z_{i+1}(\theta)$, or

,

where $i$ is the complex coefficient when multiplied and the qubit when
indexed, using $\theta=0.1$. We then do an even- and odd-layered
Trotterization for one discrete layer of the simulation. This Trotterization
on $n=5$ is

for a single discrete time-step. We use the initial condition
$\ket{\psi_{0}}=\ket{101010101010101}$.

####  D.2.2 String edit maps

The string edit maps $\mathbb{L}_{\text{HeisXXX}}$ for this simulation are
SWAP-isomorphic, with each local string edit map
$\mathcal{L}_{\text{HeisXXX}_{i,i+1}}$ having the behavior:

$\begin{split}\mathcal{L}_{\text{HeisXXX}_{i,i+1}}(\ket{00})&=\\{\ket{00}\\}\\\
\mathcal{L}_{\text{HeisXXX}_{i,i+1}}(\ket{01})&=\\{\ket{01},\ket{10}\\}\\\
\mathcal{L}_{\text{HeisXXX}_{i,i+1}}(\ket{10})&=\\{\ket{01},\ket{10}\\}\\\
\mathcal{L}_{\text{HeisXXX}_{i,i+1}}(\ket{11})&=\\{\ket{11}\\}\\\ \end{split}$
(46)

Each string edit (row in Eq. (46)) obviously conserves
$S=\sum_{j=i}^{i+1}Z_{j}$.

###  D.3 $T_{6}$ quantum cellular automata

Here we give the simulation details of the $T_{6}$ quantum cellular automata,
with model details in Sec. II.2.2 and simulation results in Sec. VI.2.

The $k$-local unitary used, ($k=3$), is
$U_{i}=P^{(1)}_{i-1}H_{i}P^{(0)}_{i+1}+P^{(0)}_{i-1}H_{i}P^{(1)}_{i+1}+P^{(1)}_{i-1}\mathds{1}_{i}P^{(1)}_{i+1}+P^{(0)}_{i-1}\mathds{1}_{i}P^{(0)}_{i+1}=\sum_{\alpha,\beta=0}^{1}P_{i-1}^{(\alpha)}\otimes(H_{i})^{\delta_{\alpha+\beta,1}}\otimes
P_{i+1}^{(\beta)}$, where $H_{i}$ is the Hadamard gate. This unitary is
implemented with the circuit

.

We implement this as a global unitary, which is applied once at each time-step
in the simulation by doing a parallel layer of all gates $U_{i}$ where $i$ is
even, then a layer of all gates where $i$ is odd.

.

We only use the initial condition $\ket{\psi_{0}}=\ket{000000010000000}$

####  D.3.1 String edit maps

The string edit maps $\mathbb{L}_{T_{6}}$ conserve the domain-wall, and each
local string edit map $\mathcal{L}_{T_{6}{i-1,i,i+1}}$ has the behavior

$\begin{split}\mathcal{L}_{T_{6}{i-1,i,i+1}}(\ket{000})&=\\{\ket{00}\\}\\\
\mathcal{L}_{T_{6}{i-1,i,i+1}}(\ket{001})&=\\{\ket{001},\ket{011}\\}\\\
\mathcal{L}_{T_{6}{i-1,i,i+1}}(\ket{010})&=\\{\ket{010}\\}\\\
\mathcal{L}_{T_{6}{i-1,i,i+1}}(\ket{011})&=\\{\ket{001},\ket{011}\\}\\\
\mathcal{L}_{T_{6}{i-1,i,i+1}}(\ket{100})&=\\{\ket{100},\ket{110}\\}\\\
\mathcal{L}_{T_{6}{i-1,i,i+1}}(\ket{101})&=\\{\ket{101}\\}\\\
\mathcal{L}_{T_{6}{i-1,i,i+1}}(\ket{110})&=\\{\ket{100},\ket{110}\\}\\\
\mathcal{L}_{T_{6}{i-1,i,i+1}}(\ket{111})&=\\{\ket{111}\\}\\\ \end{split}$
(47)

which conserves the simulation operator $S=\sum_{j=i-2}^{i+1}Z_{j}Z_{j+1}$.

###  D.4 $F_{4}$ quantum cellular automata

Here we give the simulation details of the $F_{6}$ quantum cellular automata,
with model details in Sec. II.2.3 and simulation results in Sec. VI.2 and Sec.
VI.3.

The $k$-local unitary used ($k=5$) is built on the local unitary

$U_{i}=\sum_{\alpha,\beta,\gamma,\omega=0}^{1}P_{i-2}^{\alpha}P_{i-1}^{\beta}(H_{i})^{\delta_{\alpha+\beta+\gamma+\omega,2}}P_{i+1}^{\gamma}P_{i+2}^{\omega},$
(48)

which we implement through quantum gates with the circuit

.

This local circuit is combined into a unitary applied at even discrete time-
steps

$U(\tau_{j=\text{even}})=\prod_{i=2,5,8,\ldots}U_{i}(\tau_{j})\prod_{i=3,6,9,\ldots}U_{i}(\tau_{j})\prod_{i=4,7,10,\ldots}U_{i}(\tau_{j}),$
(49)

and a unitary applied at odd discrete time-steps

$U(\tau_{j=\text{odd}})=\prod_{i=3,6,9,\ldots}U_{i}(\tau_{j})\prod_{i=2,5,8\ldots}U_{i}(\tau_{j})\prod_{i=4,7,10,\ldots}U_{i}(\tau_{j}).$
(50)

These unitaries on, for example 7 qubits, have the circuit diagrams

and

.

Notice that gates in the same product in Eq. (49) and (50), for example,
$U_{4}$ and $U_{7}$, commute because only their controls overlap. To get to,
for example $t=4$, our circuit would look like

####  D.4.1 String edit maps

Previously, no symmetry was known for the $F_{4}$ QCA. Below, we outline the
non-identity behavior of each local string edit map $\mathcal{L}_{F_{4}{i}}$:

$\begin{split}\mathcal{L}_{F_{4}{i}}(\ket{00011})&=\\{\ket{00011},\ket{00111}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{00111})&=\\{\ket{00011},\ket{00111}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{01001})&=\\{\ket{01001},\ket{01101}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{01101})&=\\{\ket{01001},\ket{01101}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{01010})&=\\{\ket{01010},\ket{01110}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{01110})&=\\{\ket{01010},\ket{01110}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{11000})&=\\{\ket{11000},\ket{11100}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{11100})&=\\{\ket{11000},\ket{11100}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{10010})&=\\{\ket{10010},\ket{10110}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{10110})&=\\{\ket{10010},\ket{10110}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{10001})&=\\{\ket{10001},\ket{10101}\\}\\\
\mathcal{L}_{F_{4}{i}}(\ket{10101})&=\\{\ket{10001},\ket{10101}\\}.\\\
\end{split}$ (51)

Any state $\ket{b}$ which does not appear above will have the behavior
$\mathcal{L}_{F_{4}i}(\ket{b})=\\{\ket{b}\\}$. It should also be understood
that $i$ is the middle qubit in the ordering above.

####  D.4.2 Initial conditions for simulations

We use a variety of initial conditions to test the imperfect path finding of
Alg. 1. Our initial condition $\ket{\psi_{0}}=\ket{000000101000000}$, tested
in Sec. VI.2, has no false rejections (in other words perfect path finding)
with $\mu=2$.

We also gather the set of symmetry-protected subspaces which have failed paths
at $\mu=9$, their average, best, and worst fidelities are shown in Sec. VI.3.
We chose an arbitrary state from each of these subspaces to be the initial
condition; the states are
$\ket{\psi_{0}}\in\\{\ket{001000011111111},\ket{011100000001000},\\\
\ket{011110101011111},\ket{011110000010001},\ket{011101000001000},\\\
\ket{011000000101010},\ket{001000011111111},\ket{011100000001000},\\\
\ket{011110101011111},\ket{011110000010001},\ket{011101000001000},\\\
\ket{011000000101010},\ket{111111101111100},\ket{110111111110101},\\\
\ket{111111101010111},\ket{111100010111110},\ket{110111111010001},\\\
\ket{110101000000100}\\}$.

## References

  * Noether [1918] E. Noether, Invariante variationsprobleme, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 1918, 235 (1918).
  * Gaillard  _et al._ [1999] M. K. Gaillard, P. D. Grannis, and F. J. Sciulli, The standard model of particle physics, Reviews of Modern Physics 71, S96 (1999).
  * Kozlov [1983] V. V. Kozlov, Integrability and non-integrability in hamiltonian mechanics, Russian Mathematical Surveys 38, 1 (1983).
  * Landau [1937] L. D. Landau, On the theory of phase transitions, Zh. Eksp. Teor. Fiz. 7, 19 (1937).
  * Fowler  _et al._ [2012] A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Surface codes: Towards practical large-scale quantum computation, Physical Review A 86, 10.1103/physreva.86.032324 (2012).
  * Bonet-Monroig  _et al._ [2018] X. Bonet-Monroig, R. Sagastizabal, M. Singh, and T. E. O’Brien, Low-cost error mitigation by symmetry verification, Phys. Rev. A 98, 062339 (2018).
  * Sagastizabal  _et al._ [2019] R. Sagastizabal, X. Bonet-Monroig, M. Singh, M. A. Rol, C. C. Bultink, X. Fu, C. H. Price, V. P. Ostroukh, N. Muthusubramanian, A. Bruno, M. Beekman, N. Haider, T. E. O’Brien, and L. DiCarlo, Experimental error mitigation via symmetry verification in a variational quantum eigensolver, Phys. Rev. A 100, 010302 (2019).
  * McClean  _et al._ [2017] J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states, Phys. Rev. A 95, 042308 (2017).
  * Cai [2021] Z. Cai, Multi-exponential error extrapolation and combining error mitigation techniques for NISQ applications, npj Quantum Information 7, 10.1038/s41534-021-00404-3 (2021).
  * Larocca  _et al._ [2022] M. Larocca, F. Sauvage, F. M. Sbahi, G. Verdon, P. J. Coles, and M. Cerezo, Group-invariant quantum machine learning, PRX Quantum 3, 030341 (2022).
  * Daley  _et al._ [2022] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Practical quantum advantage in quantum simulation, Nature 607, 667 (2022).
  * Acharya  _et al._ [2022] R. Acharya, I. Aleiner, R. Allen, T. I. Andersen, M. Ansmann, F. Arute, K. Arya, A. Asfaw, J. Atalaya, R. Babbush, _et al._ , Suppressing quantum errors by scaling a surface code logical qubit, arXiv preprint arXiv:2207.06431 (2022).
  * Cai  _et al._ [2022] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O’Brien, Quantum error mitigation, arXiv preprint arXiv:2210.00921 (2022).
  * Quek  _et al._ [2022] Y. Quek, D. S. França, S. Khatri, J. J. Meyer, and J. Eisert, Exponentially tighter bounds on limitations of quantum error mitigation, arXiv preprint arXiv:2210.11505 (2022).
  * Takagi  _et al._ [2022] R. Takagi, S. Endo, S. Minagawa, and M. Gu, Fundamental limits of quantum error mitigation, npj Quantum Information 8, 10.1038/s41534-022-00618-z (2022).
  * Arute  _et al._ [2020] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, A. Bengtsson, S. Boixo, M. Broughton, B. B. Buckley, _et al._ , Observation of separated dynamics of charge and spin in the fermi-hubbard model, arXiv preprint arXiv:2010.07965 (2020).
  * Hébert  _et al._ [2001] F. Hébert, G. G. Batrouni, R. T. Scalettar, G. Schmid, M. Troyer, and A. Dorneich, Quantum phase transitions in the two-dimensional hardcore boson model, Physical Review B 65, 014513 (2001).
  * Anschuetz  _et al._ [2022] E. R. Anschuetz, A. Bauer, B. T. Kiani, and S. Lloyd, Efficient classical algorithms for simulating symmetric quantum systems, arXiv preprint arXiv:2211.16998 (2022).
  * Radjavi and Rosenthal [2003] H. Radjavi and P. Rosenthal,  _Invariant subspaces_ (Courier Corporation, 2003).
  * Suzuki [1991] M. Suzuki, General theory of fractal path integrals with applications to many-body theories and statistical physics, Journal of Mathematical Physics 32, 400 (1991).
  * Jones  _et al._ [2022] E. B. Jones, L. E. Hillberry, M. T. Jones, M. Fasihi, P. Roushan, Z. Jiang, A. Ho, C. Neill, E. Ostby, P. Graf, E. Kapit, and L. D. Carr, Small-world complex network generation on a digital quantum processor, Nature Communications 13, 4483 (2022).
  * Hillberry  _et al._ [2021] L. E. Hillberry, M. T. Jones, D. L. Vargas, P. Rall, N. Y. Halpern, N. Bao, S. Notarnicola, S. Montangero, and L. D. Carr, Entangled quantum cellular automata, physical complexity, and goldilocks rules, Quantum Science and Technology 6, 045017 (2021).
  * Lidl [1998] R. Lidl,  _Applied abstract algebra_ , second edition. ed., Undergraduate Texts in Mathematics (Springer, New York, 1998) p. 337.
  * Fischer and Meyer [1971] M. J. Fischer and A. R. Meyer, Boolean matrix multiplication and transitive closure, in  _12th Annual Symposium on Switching and Automata Theory (swat 1971)_ (1971) pp. 129–131.
  * Skiena [2008] S. S. Skiena, Sorting and searching, in  _The Algorithm Design Manual_ (Springer London, London, 2008) p. 477.
  * Lee [1961] C. Y. Lee, An algorithm for path connections and its applications, IRE Transactions on Electronic Computers EC-10, 346 (1961).
  * Kullback and Leibler [1951] S. Kullback and R. A. Leibler, On information and sufficiency, The annals of mathematical statistics 22, 79 (1951).
  * Kapit  _et al._ [2020] E. Kapit, P. Roushan, C. Neill, S. Boixo, and V. Smelyanskiy, Entanglement and complexity of interacting qubits subject to asymmetric noise, Phys. Rev. Res. 2, 043042 (2020).
  * Endo  _et al._ [2021] S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, Hybrid quantum-classical algorithms and quantum error mitigation, Journal of the Physical Society of Japan 90, 032001 (2021).
  * Needleman and Wunsch [1970] S. B. Needleman and C. D. Wunsch, A general method applicable to the search for similarities in the amino acid sequence of two proteins, Journal of Molecular Biology 48, 443 (1970).