# Matchgate benchmarking: Scalable benchmarking of a continuous family of
many-qubit gates

Jonas Helsen  QuSoft & Korteweg-de Vries Institute for Mathematics, University
of Amsterdam, Science Park 123, 1098 XG Amsterdam, The Netherlands  Centrum
Wiskunde & Informatica (CWI), Science Park 123, 1098 XG Amsterdam, The
Netherlands  Sepehr Nezami  Institute for Quantum Information and Matter,
Caltech, Pasadena, CA 91125, USA  Matthew Reagor  Rigetti Computing, 775 Heinz
Ave, Berkeley, CA 94710, USA  Michael Walter  QuSoft & Korteweg-de Vries
Institute for Mathematics, University of Amsterdam, Science Park 123, 1098 XG
Amsterdam, The Netherlands  Institute for Theoretical Physics & ILLC,
University of Amsterdam, Science Park 123, 1098 XG Amsterdam, The Netherlands
Faculty of Computer Science, Ruhr University Bochum, Universitätsstraße 150,
44801 Bochum, Germany

###### Abstract

We propose a method to reliably and efficiently extract the fidelity of many-
qubit quantum circuits composed of continuously parametrized two-qubit gates
called matchgates. This method, which we call _matchgate benchmarking_ ,
relies on advanced techniques from randomized benchmarking as well as insights
from the representation theory of matchgate circuits. We argue the formal
correctness and scalability of the protocol, and moreover deploy it to
estimate the performance of matchgate circuits generated by two-qubit XY spin
interactions on a quantum processor.

Quantum computers promise a revolution in computational power, and a
multinational effort is underway to construct them. One of the key challenges
in the building and operating of quantum computers is the appearance of errors
in computations, either due to inaccuracies in control or due to interactions
with the environment. It is thus vitally important to be able to characterize
accurately and efficiently the type and magnitude of errors present in quantum
operations. To this end, a variety of techniques have been developed, with the
most popular class of techniques known as randomized benchmarking (RB) [1, 2,
3, 4, 5, 6, 7, 8, 9], where one characterizes the quality of gates in a
gateset by applying random sequences of gates of increasing length, and tracks
the corresponding increase in average error. For a recent overview of RB
protocols see [10] and references therein.

Randomized benchmarking has been extremely successful in characterizing
quantum operations on a variety of platforms [11, 6, 12]. Yet it suffers from
a number of shortcomings that limit its usefulness in some important
situations. Firstly, standard RB protocols mix the error contributions of
various types of gates (such as single qubit gates and two qubit gates) and
only report an average error. This is problematic because different types of
gates are created by different physical mechanisms and hence have different
error contributions. Moreover, different types of gates contribute differently
to the error thresholds that must be met for fault-tolerance. For instance it
is often the case that more stringent requirements are imposed on two-qubit
gates than on single-qubit gates. Secondly, standard randomized benchmarking
protocols only test discrete gatesets (such as the Clifford group), while
continuously parametrized gatesets are vital for near-term quantum computing
applications such as VQE and QAOA [13, 14]. For these reasons it is desirable
to devise gate assessment procedures that combine the proven advantages of
randomized benchmarking with the ability to handle continuous gate families
and focus on a single type of quantum gates.

In this paper, we address this challenge by proposing _matchgate benchmarking_
, an advanced randomized benchmarking procedure based on the general framework
given in [10] as well as the recently introduced linear cross-entropy
benchmarking [15]. Our procedure natively uses continuously parametrized two-
qubit gates and estimates fidelities in a scalable way, both in terms of
statistical sampling and classical computational resources. To prove its
correctness and scalability we use techniques from the representation theory
of the matchgate group. Moreover we provide an implementation of the protocol
on a small quantum computer, showing that our protocol can reliably assess the
quality of two qubit gates in a realistic environment.

Because this paper presents both proposals for experimental procedures and
mathematical results, and thus is aimed at both experimental practitioners and
theorists, we defer all technical proofs to appendices.

##  1 Matchgates

Matchgate circuits are a continuous class of quantum circuits originally
conceived by Valiant [16] and subsequently connected to the theory of free
fermions by Knill [17] and Terhal-DiVincenzo [18], see also [19, 20, 21]. They
are explicitly realized by $\operatorname{XY}$ or by $\operatorname{XX}$ and
$\operatorname{YY}$ spin interactions and are thus the natural choice for two-
qubit gates on many physical platforms such as ion traps. The standard
$\operatorname{iSWAP}$ [22] and $\operatorname{XY}(\theta)$ gates [23] are
examples of matchgates (though they do not generate the full group by
themselves). A key property of matchgate circuits is that they are efficiently
classically simulable (like the better known Clifford group often used in
standard RB), which is a necessary requirement for a scalable randomized
benchmarking procedure. The two-qubit matchgates are generated by unitaries of
the form $U(\alpha)=\exp(i\alpha P\otimes P^{\prime})$, where $P,P^{\prime}$
are Pauli $X$ or $Y$ matrices. The $n$-qubit _matchgate group_
$\mathcal{M}_{n}$ is then defined by considering $n$ qubits on a line and
composing _nearest-neighbor_ gates of this form, along with single-qubit
$Z(\theta)=\exp(i\theta Z)$ gates. We further extend this group with a Pauli
$X$ on the _last_ qubit. This forms the
$\mathcal{M}^{+}_{n}=\langle\mathcal{M}_{n},X_{n}\rangle$ group, which we will
refer to as _generalized matchgates_.

Matchgates are intimately connected to non-interacting fermions. This
connection is key to their efficient simulation. To see this, consider the
Majorana fermion operators (represented on qubits by the Jordan-Wigner
isomorphism)

$\displaystyle\gamma_{2j-1}$ $\displaystyle=Z_{1}\ldots
Z_{j-1}X_{j}I_{j+1}\ldots I_{n},$ $\displaystyle\gamma_{2j}$
$\displaystyle=Z_{1}\ldots Z_{j-1}Y_{j}I_{j+1}\ldots I_{n},$

with $j\in[n]:=\\{1,\dots,n\\}$ and $I,X,Y,Z$ the Pauli matrices, the
subscript indicating on which qubit they act. Any matchgate unitary
$U\in\mathcal{M}_{n}$ acts by conjugation on the Majorana operators as

$U\gamma_{j}U^{\dagger}=\sum_{j\in[2n]}R_{kj}\gamma_{k},$ (1)

with $R\in\operatorname{SO}(2n)$ a rotation matrix. Moreover the $X_{n}$ gate
maps all $\gamma_{i}$ for $i<2n$ to themselves but maps $\gamma_{2n}$ to
$-\gamma_{2n}$, so it corresponds to a reflection $F$ of the $2n$-th axis. In
this way, matchgate and generalized matchgates unitaries can be efficiently
tracked on a classical computer. Moreover, for any rotation there is a
matchgate unitary $U=U(R)$ implementing it, and for every element $Q$ of
$\operatorname{O}(2n)$ (rotations plus reflections) there is a corresponding
generalized matchgate unitary $U=U(Q)$. In fact, $\mathcal{M}_{n}$ is
generated as a Lie group by Hamiltonians of the form
$H=\frac{i}{4}\sum_{j,k\in[2n],j\neq k}\alpha_{jk}\gamma_{j}\gamma_{k}$, with
$\alpha$ a real antisymmetric $2n\times 2n$ matrix, so $\mathcal{M}_{n}$ can
be understood as a representation of $\operatorname{Spin}(2n)$ and
$\mathcal{M}_{n}^{+}$ a representation of $\operatorname{Pin}(2n)$. Since the
Majorana operators and their products span the space of $n$-qubit operators,
(generalized) matchgate unitaries are fully determined by the corresponding
rotation (and reflection) matrix (up to an overall phase).

Before we can define our benchmarking procedure we need to briefly discuss the
action of the generalized matchgate group on the space of $n$-qubit operators.
We denote products of Majoranas as $\gamma[S]=\prod_{s\in S}\gamma_{s}$ for
$S\subseteq[2n]$, with $\gamma_{\emptyset}=I$ and the product taken in
increasing order. For each $k\in[2n]$, consider the subspace
$\Gamma_{k}=\langle\gamma[S]\;|\;S\subseteq[2n],\;\lvert S\rvert=k\rangle$
spanned by products of $k$ Majorana operators. Then, for $k\in[2n]$ each
$\Gamma_{k}$ is an irreducible representation of the generalized matchgate
group $\mathcal{M}_{n}^{+}$. Moreover all these representations are
inequivalent. A proof of this statement is given in Lemma 3.

We note that the addition of the extra bit flip gate, lifting the matchgates
to the generalized matchgates, is critical in ensuring the mutual
inequivalence of the representations $\Gamma_{k}$.

##  2 Matchgate benchmarking

The matchgate benchmarking protocol, given formally in Algorithm 1, estimates
the quality of generic circuits in the generalized matchgate group
$\mathcal{M}_{n}^{+}$ in a manner that scales efficiently with the number of
qubits and is resistant to state preparation and measurement (SPAM) errors.
The output of matchgate benchmarking is a list of decay parameters
$\lambda_{k}$ that characterize the noise associated to the subspace
$\Gamma_{k}$. We will call these decay parameters _Majorana fidelities_.

The protocol consists of multiple rounds with varying parameters $k\in[2n]$
and sequence lengths $m$. Each round starts with the preparation of either the
all-zero $\lvert 0_{n}\rangle:=\lvert 0\rangle^{\otimes n}$ or the all-plus
$\lvert+_{n}\rangle:=\lvert+\rangle^{\otimes n}$ state. This is followed by
$m$ generalized matchgate unitaries $U(Q_{1}),\dots,U(Q_{m})$, chosen
uniformly and independently at random from $\mathcal{M}_{n}^{+}$ (we describe
an efficient method for sampling below). Finally, all qubits are measured in
either the computational ($Z$) basis or the Hadamard ($X$) basis. We write
$\rho_{0}$ for the initial state and $\\{E_{x}\\}_{x\in\\{0,1\\}^{n}}$ for the
measurement POVM, and refer to the two SPAM settings as $X$ and $Z$-basis
SPAM. The above is repeated many times until the relative frequencies $f_{x}$
of the measurement outcomes $x\in\\{0,1\\}^{n}$ give a good estimate of the
true probabilities, which we denote by $p(x|Q,m)$, where $Q=Q_{m}\cdots
Q_{1}$. Since $Q$ itself is uniformly random, by averaging over many random
sequences we thus obtain a good estimate $\hat{f}_{k}(m)$ of the weighted
average

$\displaystyle
f_{k}(m)=\int_{\operatorname{O}(2n)}dQ\sum_{x\in\\{0,1\\}^{n}}\,\alpha_{k}(x,Q)p(x|Q,m),$
(2)

where we use the _correlation function_ $\alpha_{k}$ defined by

$\displaystyle\alpha_{k}(x,Q)=\frac{1}{N_{k}}\\!\operatorname{Tr}\Bigl{(}E_{x}P_{k}\bigl{(}U(Q)\rho_{0}U(Q)^{\dagger}\bigr{)}\Bigr{)}.$
(3)

Here $P_{k}$ denotes the projection superoperator onto the subspace
$\Gamma_{k}$, and the normalization constant $N_{k}=2^{-n}\binom{n}{\lfloor
k/2\rfloor}^{2}\binom{2n}{k}^{-1}$ is chosen so that $f_{k}(m)=1$ if the gates
are perfectly implemented. The correlation functions $\alpha_{k}$ can be
efficiently computed using the classical simulation techniques for matchgates
[18, 17] as well as several tricks for evaluating Pfaffian sums [24]. We give
explicit expressions in Appendix E.

##  3 Interpretation and analysis

Intuitively the randomization over gates $U(Q_{1}),\ldots,U(Q_{m})$ averages
out the noise associated to each gate, resulting in a linear combination of
generalized depolarizing channels, one for each irreducible subspace
$\Gamma_{k}$. This observation forms the basis of the randomized benchmarking
approach, of which matchgate benchmarking is an example 111For comparison: In
standard RB with the Clifford group there are two subrepresentations, giving
rise to a standard depolarizing channel upon averaging.. In our setting, the
associated depolarization parameters are described by what we call the
_Majorana fidelities_ $\lambda_{k}$, since the corresponding subspaces
$\Gamma_{k}$ are spanned precisely by the $k$-fold product of the Majorana
operators. The role of the correlation functions $\alpha_{k}$ is precisely to
address the individual subspaces $\Gamma_{k}$ and thus isolate the
corresponding Majorana fidelities. In particular $f_{k}(m)$ will be determined
precisely by the $m$-th power $\lambda_{k}^{m}$. To make this more concrete,
suppose that each generalized matchgate unitary $U(Q)$ is realized by a
quantum channel $\Phi(Q)$ that describes its actual implementation.222The
existence of such a map $\Phi$ is an assumption on the underlying device and
in particular excludes time- and context-dependent effects. It is, however,
the weakest assumption under which RB protocols can be guaranteed to function
correctly [10]. If we assume _gate-independent noise_ , i.e., that there is a
quantum channel $\Lambda$ such that $\Phi(Q)(\rho)=\Lambda(U(Q)\rho
U(Q)^{\dagger})$ for all $Q\in\operatorname{O}(2n)$, then the average
$f_{k}(m)$ is described _exactly_ by a single exponential decay

$f_{k}(m)=A_{k}\lambda_{k}^{m}.$

Importantly, $\lambda_{k}$ depends only on the noise channel $\Lambda$, while
$A_{k}$ also depends on the SPAM. A proof of of the above statement is given
in Appendix C The assumption of gate-independent noise is unrealistic, but it
can be relaxed significantly using Fourier analytic techniques [10]. We
outline an approach to this in Appendix D but leave a detailed derivation for
future work.

In the absence of noise, the Majorana fidelities are equal to one, so their
deviation from the identity encodes properties of the noise. In particular,
the average gate fidelity can (for gate-independent noise) be recovered by

$2^{-n}\sum_{k=0}^{2n}\\!\binom{2n}{k}\lambda_{k}=\\!(2^{n}\\!+\\!1)F_{\mathrm{avg}}(\Lambda)-1.$
(4)

Interpreting the individual values $\lambda_{k}$ operationally is less
straightforward. The parameter $\lambda_{0}$ has a well-known interpretation
as a measure of trace-loss of the channel $\Lambda$. Moreover, if $\Lambda$ is
unitary (i.e. generated by some Hamiltonian $H$) then the parameter
$\lambda_{2n}$ can be seen as a measure of parity preservation of this
evolution. More generally we can interpret the parameters $\lambda_{k}$ if we
restrict the channel $\Lambda$ to be Gaussian (as discussed in [20]) and
unital. In this case the channel $\Lambda$ has associated to it a $2n\times
2n$ real matrix $B$ s.t. $BB^{T}\leq\mathbb{1}$ and the action on a Majorana
operator $\gamma[S]$ is defined as
$\Lambda(\gamma[S])=\sum_{S^{\prime}\subset[2n]\,,|S^{\prime}|=|S|}\det(B[S^{\prime},S])\gamma[S^{\prime}]$,
where $B[S^{\prime},S]$ denotes the submatrix of $B$ with rows in $S^{\prime}$
and columns in $S$. Hence the decay parameter $\lambda_{k}$ is in this case
precisely given by the average over the principal minors of size $k$. It is
easy to see that these averages over principal minors precisely correspond to
(normalized) elementary symmetric polynomials of the eigenvalues of $B$. Hence
the Majorana fidelities encode eigenvalue information of Gaussian noise. These
eigenvalues can then be in principle extracted by solving the system of $2n$
polynomials in $2n$ unknowns (although one would have to work out the
stability of solutions under statistical noise, which we do not attempt here).

Given the functional form of $f_{k}(m)$, we can extract the value of
$\lambda_{k}$ by fitting to a single exponential decay. To perform this fit in
practice, a correct choice of SPAM operators $\rho_{0},\\{E_{x}\\}$ is
required that ensures that the pre-factors $A_{k}$ are large enough. In the
noise-free limit the parameters $A_{k}$ can be explicitly computed (which we
do in Appendix C). For odd $k$ (and $X$-basis SPAM), we find that $A_{k}=1$,
and thus the fitting problem is well conditioned. For even $k$ (and $Z$-basis
SPAM), we similarly find that $A_{k}=1$. We also note that if we choose
$Z$-basis SPAM and consider odd $k$, we have $A_{k}=0$, giving no visibility.
Similarly we have $A_{2n}=0$ for $X$-basis SPAM. This motivates our use of
different SPAM settings for different values of $k$.

Algorithm 1 Matchgate benchmarking

1:for $k\in\\{0,\dots,2n\\}$ do

2:for $m\in\text{sequence lengths}$ do

3:for $i\in[K]$ do

4: Prepare $\lvert 0\rangle^{\otimes n}$ ($k$ even) or
$\lvert+\rangle^{\otimes n}$ ($k$ odd).

5: for $j\in[m]$ do

6: Apply $Q_{j}^{(i)}\in\operatorname{O}(2n)$ chosen uniformly at random.

7: end for.

8: Measure in the $Z$ ($k$ even) or $X$ basis ($k$ odd).

9: Repeat the above many times and record frequencies $f^{(i)}_{x}$ of
measurement outcomes $x\in\\{0,1\\}^{n}$.

10:end for

11:Compute the empirical average
$\displaystyle\\!\\!\\!\\!\hat{f}_{k}(m)=\frac{1}{K}\sum_{i=1}^{K}\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q_{m}^{(i)}\cdots
Q_{1}^{(i)})f^{(i)}_{x}.$ (5)

12:end for

13:Fit $\\{\hat{f}_{k}(m)\\}_{m}$ to
$\hat{f}_{k}(m)=_{\mathrm{fit}}A_{k}\lambda^{m}_{k}$.

14:end for

15:Output the Majorana fidelities $\\{\lambda_{k}\\}_{k=0}^{2n}$.

##  4 Generating random matchgate circuits

Random rotations in $\operatorname{O}(2n)$ (and thus generalized matchgate
unitaries) can readily be sampled in an efficient manner by a variety of
methods. However it is desirable to generate them directly in the form of
circuits involving one- and two-qubit matchgates. Begin by noting that we can
decompose any generalized matchgate $U(Q)$ as $U(Q)=U(R)X_{n}^{b}$ with
$b\in\\{0,1\\}$ where $U(R)$ (with $R\in\operatorname{SO}(2n)$) is a matchgate
and $X_{n}$ is a bit-flip on the last qubit. Hence the task of sampling random
generalized matchgates reduces to that of sampling random matchgates. For this
we give a method based on the probabilistic Hurwitz lemma [25]. Consider the
rotation

$\displaystyle
R=(G^{(1)}_{2n\\!-\\!1}\\!\cdots\\!G^{(1)}_{1})(G^{(2)}_{2n\\!-\\!1}\\!\cdots\\!G^{(2)}_{2})\\!\cdots\\!G^{(2n\\!-\\!1)}_{2n\\!-\\!1},$
(6)

where each $G^{(i)}_{j}$ is a two-dimensional rotation by a random angle in
the $j,j+1$-plane. Proposition 2.1 in [25] implies that $R$ is a uniformly
(Haar) random rotation in $\operatorname{SO}(2n)$. The formula for $R$
translates directly into a circuit for the corresponding matchgate unitary
$U(R)$ that only involves single qubit $Z_{j}$ and two-qubit $X_{j}X_{j+1}$
rotations. Since $Z$-rotations are virtual [26, 11], and hence noiseless, in
many platforms, the dominant source of noise in this circuit is from the two
qubit gates. Adding the aforementioned random $X_{n}$ gates we obtain a
uniformly random circuit $U(Q)\in\mathcal{M}_{n}^{+}$ generalized matchgate
circuit $U(Q)$ given above is comparable to that of a generic $n$-qubit
Clifford gate [27].

Finally, if one only has access to $\operatorname{XY}$ gates (corresponding to
$\operatorname{XX}+\operatorname{YY}$ rotations, which do not by themselves
generate the full matchgate group) as opposed to $\operatorname{XX}$ or
$\operatorname{YY}$ rotations, the above construction can be implemented by
the identity
$\operatorname{XX}(\theta)=\operatorname{XY}(\theta/2)\,X_{1}\,\operatorname{XY}(\theta/2)\,X_{1}$,
where $X_{1}$ denotes the bit-flip gate on the first qubit.

We conjecture that one can also efficiently sample _approximately uniform_
matchgate unitaries by repeatedly choosing nearest-neighbor pairs of qubits at
random and applying a random element of $\mathcal{M}_{2}$. This is a variation
of the well-known Kac random walk on $\operatorname{SO}(2n)$ which mixes to
approximate uniformity in polynomial time [28]. This maybe possibly be an even
more efficient way of sampling (approximately) uniformly random matchgate
circuits.

##  5 Statistical scalability

We now consider the scalability of the matchgate benchmarking protocol with
respect to the number of qubits. Recall that in Algorithm 1 we determine the
relative frequencies of measurement outcomes for a number of random matchgate
sequences. By taking an empirical average, one obtains an estimate
$\hat{f}_{k}(m)$ for $f_{k}(m)$. It is a priori unclear whether the variance
of this estimate might grow exponentially with the number of qubits $n$,
rendering the estimation procedure infeasible beyond a few qubits. We argue
this is not the case by explicitly bounding the variance in the noise-free
limit.

######  Theorem 1.

Consider the estimator $\hat{f}_{k}(m)$ for the quantity $f_{k}(m)$ defined in
Eqs. 2 and 5. Assuming no noise, its variance is bounded (uniformly in $k$ and
$m$) as

$\mathds{V}(\hat{f}_{k}(m))=\frac{1}{K}O(\mathrm{poly}(n)).$

A proof of this theorem is given in Appendix A. The central ingredient in this
theorem is a novel moment bound for random matchgates based on the
representation theory of $\operatorname{SO}(2n)$, which may be of independent
interest:

######  Lemma 2.

Let $\lvert\theta\rangle$ be the all-zero ($\lvert 0\rangle^{\otimes n}$) or
the all-plus ($\lvert+\rangle^{\otimes n}$) state, and let $t$ be a fixed
integer. Then:

$\int_{\operatorname{SO}(2n)}|\langle\theta\rvert
U(Q)\lvert\theta\rangle|^{2t}=2^{-tn}\,O(\mathrm{poly}(n)).$

A proof of this lemma is given in Appendix B. We note that our variance upper
bound in Theorem 1 is likely loose and we expect the real variance to be
substantially smaller. The theorem can also be extended to the case of gate-
dependent noise at the cost of some technical complications, but we do not
pursue this here.

##  6 Experimental demonstration

We apply the matchgate benchmarking protocol (Algorithm 1) to benchmark the
native $XY(\theta)$ gate [23] between two qubits on the Rigetti Aspen-8 chip.
A collection of circuits was run consecutively with fixed sequence lengths
ranging from $m=2$ to $m=24$, sampling $K=64$ random sequences of orthogonal
matchgates per sequence length, and performing $L=400$ measurement repetitions
(shots) per sequence; for a total of $5.9\times 10^{5}$ individual shots. The
code and data for this experiment can be found at [29]. Figure 1 shows the
results of this experiment. All error bars are bootstrapped $95\%$ confidence
intervals. For $n=2$, there are five exponential decays, associated to
$k\in[4]$. From fitting the experimental data to single exponentials we obtain
the following values for the Majorana fidelities:

$\lambda_{0}$ | ${1.000\\!\pm\\!0.001}$
---|---
$\lambda_{1}$ | ${0.78\\!\pm\\!0.05}$
$\lambda_{2}$ | ${0.85\\!\pm\\!0.02}$
$\lambda_{3}$ | ${0.87\\!\pm\\!0.02}$
$\lambda_{4}$ | ${0.83\\!\pm\\!0.02}$

From these decays and Eq. 4 we can infer that the expected average fidelity of
a random two-qubit matchgate is $F={0.88\pm 0.02}$. Since a random two qubit
matchgate requires four $XY(\theta)$ gates, we can make a heuristic lower
bound estimate of the average fidelity of the $XY(\theta)$ gate by assuming
that the single-qubit gates are essentially noiseless and that the fidelity is
approximately multiplicative, finding $F={0.97}$ for a single $XY(\theta)$
gate. We compare this estimate to standard interleaved RB applied to the same
$\operatorname{iSWAP}$ gate ($XY(\pi)$), where we observe $F={0.986\pm
0.006}$, consistent with the fidelity range for $XY$-gates observed in [23].
We attribute additional error in the matchgate construction to single-qubit
rotations.

Figure 1: Five exponential decays associated with performing matchgate
benchmarking using circuits generated by $XY(\theta)$ gates. For readability
each exponential decay is offset vertically by some amount. Based on this data
we can conclude a two-qubit matchgate circuit fidelity of
$F_{\mathrm{avg}}={0.88\pm 0.02}$.

##  7 Extensions and variations of the gateset

The benchmarking procedure given in Algorithm 1 for nearest-neighbor
generalized matchgate circuits on a line can readily be adapted to related
gatesets with interesting properties.

First, following Knill [17] (see also [30]), one can extend the generalized
matchgates by arbitrary single-qubit gates on the first qubit on the line.
This corresponds to gates generated by Hamiltonians composed of linear
Majorana terms, i.e. $H_{\mathrm{lin}}=\sum_{i}v_{i}\gamma_{i}$. Note that the
distinction between matchgates and generalized matchgates collapses in this
case. Equivalently, one can add rotations along the $ZX$ axis between the
first two neighboring qubits (a cross-resonance gate [31, 32]). This extension
corresponds to the group $\operatorname{SO}(2n+1)$ [17], which has the spaces
$\Gamma_{0}$ and $\Gamma_{2k^{\prime}-1}\oplus\Gamma_{2k^{\prime}}$ for
$k^{\prime}\in[n]$ as mutually inequivalent irreducible subspaces. Matchgate
benchmarking is easily adapted by using the correlation functions
$\alpha^{\prime}_{k^{\prime}}(x,R)\propto\operatorname{Tr}\bigl{(}E_{x}(P_{2{k^{\prime}}-1}\\!+\\!P_{2k^{\prime}})(U(R)\rho_{0}U(R)^{\dagger})\bigr{)}$
with $R\in\operatorname{SO}(2n+1)$. Assuming gate-independent noise, one finds
that $f_{k^{\prime}}(m)=A_{k^{\prime}}\lambda_{k^{\prime}}^{m}$, where the
$n+1$ parameters $\lambda_{k^{\prime}}\in\mathbb{R}$ depend only on the noisy
implementation of the circuits, and the average gate fidelity can be extracted
as
$2^{-n}\sum_{k^{\prime}=0}^{n}\binom{2n+1}{k^{\prime}}\lambda_{k^{\prime}}=(2^{n}+1)F_{\mathrm{avg}}(\Lambda)-1$.

Second, we can also extend the nearest-neighbor generalized matchgate circuits
on a line to those on a circle. This corresponds to gates generated by
Hamiltonians of the form
$H_{\mathrm{circle}}=\sum_{i,j}\alpha_{i,j}\gamma[\\{i,j\\}]+\sum_{i,j}\beta_{i,j}\gamma[\\{i,j\\}^{c}]$,
together with a single $X$ gate. Their classical simulability was, to our
knowledge, first noted in [33]. One can again work out the corresponding
representation theory and write down appropriate correlation functions.

Third, one can also perform matchgate benchmarking with the ordinary matchgate
group $\mathcal{M}_{n}$ (without the additional bit-flip gate). Now the
representations $\Gamma_{k},\Gamma_{2n-k}$ for $k<n$ become equivalent, and
the representation $\Gamma_{n}$ splits into two inequivalent representations.
The correlation functions change to
$\alpha^{\prime\prime}_{k}(x,R)\propto\operatorname{Tr}\bigl{(}E_{x}(P_{k}\\!+\\!P_{2n-k})(U(R)\rho_{0}U(R)^{\dagger})\bigr{)}$.
Due to the equivalence of representations, the data $f_{k}(m)$ for $k\in[n]$
must be fitted to a $2\times 2$ _matrix_ -exponential decay
$\operatorname{Tr}(A_{k}M_{k}^{m})$, with the eigenvalues of the matrices
$M_{k}$ carrying the fidelity information. This is a significantly harder
fitting problem in practice.

Finally, we note that the (orthogonal) matchgate group can be conjugated by a
Clifford operator (i.e. consider gates of the form $CU(Q)C^{\dagger}$ where
$U(Q)\in\mathcal{M}_{n}$ and $C$ a Clifford operator) without losing classical
simulability [21]. This conjugation leaves the representation structure, and
hence the matchgate benchmarking protocol, unchanged (apart from an
appropriate change of initial states and measurement basis). One can for
instance consider the orthogonal matchgate group rotated by single qubit
Hadamard gates (on each qubit). This rotated orthogonal matchgate group is
generated by nearest neighbor $ZZ$ rotations, single qubit $X$ rotation, and
$Z$ phase flips. As the $ZZ$ interaction is natural in superconducting circuit
qubits and often used to generate two qubit gates [34, 35], this potentially
extends the usefulness of matchgate benchmarking.

We would like to acknowledge Harold Nieuwboer, Sergii Strelchuk, Ingo Roth,
and Emilio Onorati for useful conversations. MW acknowledges support by an NWO
Veni Innovational Research Grant no. 680-47-459, NWO grant OCENW.KLEIN.267,
and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
under Germany’s Excellence Strategy - EXC 2092 CASA - 390781972. This material
is based upon work supported by the Defense Advanced Research Projects Agency
(DARPA) under agreement No. HR00112090058. While preparing this manuscript, we
became aware of [36], where a similar procedure for the standard matchgate
group is proposed.

## References

  * Hashagen et al. [2018] A. K. Hashagen, S. T. Flammia, D. Gross, and J. J. Wallman.  Real randomized benchmarking.  _Quantum_ , 2:85, 2018.  doi: 10.22331/q-2018-08-22-85.
  * Helsen et al. [2019] J. Helsen, X. Xue, L. M. K. Vandersypen, and S. Wehner.  A new class of efficient randomized benchmarking protocols.  _npj Quant. Inf._ , 5:1–9, 2019.  doi: 10.1038/s41534-019-0182-7.
  * Cross et al. [2016] A. W. Cross, E. Magesan, L. S. Bishop, J. A. Smolin, and J. M. Gambetta.  Scalable randomised benchmarking of non-Clifford gates.  _npj Quant. Inf._ , 2:16012, 2016.  doi: 10.1038/npjqi.2016.12.
  * Carignan-Dugas et al. [2015] A. Carignan-Dugas, J. J. Wallman, and J. Emerson.  Characterizing universal gate sets via dihedral benchmarking.  _Phys. Rev. A_ , 92:060302, 2015.  doi: 10.1103/PhysRevA.92.060302.
  * Wallman et al. [2015] J. J. Wallman, M. Barnhill, and J. Emerson.  Robust characterization of loss rates.  _Phys. Rev. Lett._ , 115:060501, 2015.  doi: 10.1103/PhysRevLett.115.060501.
  * Barends et al. [2014] R. Barends, J. Kelly, A. Veitia, A. Megrant, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, I.-C. Hoi, E. Jeffrey, C. Neill, P. J. J. O’Malley, J. Mutus, C. Quintana, P. Roushan, D. Sank, J. Wenner, T. C. White, A. N. Korotkov, A. N. Cleland, and John M. Martinis.  Rolling quantum dice with a superconducting qubit.  _Phys. Rev. A_ , 90:030303, 2014.  doi: 10.1103/PhysRevA.90.030303.
  * Gambetta et al. [2012] J. M. Gambetta, A. D. Córcoles, S. T. Merkel, B. R. Johnson, J. A. Smolin, J. M. Chow, C. A. Ryan, C. Rigetti, S. Poletto, T. A. Ohki, M. B. Ketchen, and M. Steffen.  Characterization of addressability by simultaneous randomized benchmarking.  _Phys. Rev. Lett._ , 109:240504, 2012.  doi: 10.1103/PhysRevLett.109.240504.
  * Knill et al. [2008] E. Knill, D. Leibfried, R. Reichle, J. Britton, R. B. Blakestad, J. D. Jost, C. Langer, R. Ozeri, S. Seidelin, and D. J. Wineland.  Randomized benchmarking of quantum gates.  _Phys. Rev. A_ , 77:012307, 2008.  doi: 10.1103/PhysRevA.77.012307.
  * França et al. [2021] Daniel Stilck França, Sergii Strelchuk, and Michał Studziński.  Efficient classical simulation and benchmarking of quantum processes in the Weyl basis.  _Physical Review Letters_ , 126(21):210502, 2021\.  doi: 10.1103/PhysRevLett.126.210502.
  * Helsen et al. [2020] Jonas Helsen, Ingo Roth, Emilio Onorati, Albert Werner, and Jens Eisert.  A general framework for randomized benchmarking.  arXiv:2010.07974, 2020.
  * Xue et al. [2019] X Xue, TF Watson, J Helsen, Daniel R Ward, Donald E Savage, Max G Lagally, Susan N Coppersmith, MA Eriksson, S Wehner, and LMK Vandersypen.  Benchmarking gate fidelities in a Si/SiGe two-qubit device.  _Phys. Rev. X_ , 9(2):021011, 2019.  doi: 10.1103/PhysRevX.9.021011.
  * Erhard et al. [2019] A. Erhard, J. J. Wallman, L. Postler, M. Meth, R. Stricker, E. A. Martinez, P. Schindler, T. Monz, J. Emerson, and R. Blatt.  Characterizing large-scale quantum computers via cycle benchmarking.  _Nature Comm._ , 10, 2019.  doi: 10.1038/s41467-019-13068-7.
  * McClean et al. [2016] Jarrod R McClean, Jonathan Romero, Ryan Babbush, and Alán Aspuru-Guzik.  The theory of variational hybrid quantum-classical algorithms.  _New J. Phys._ , 18(2):023023, 2016.  doi: 10.1088/1367-2630/18/2/023023.
  * Farhi et al. [2014] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  arXiv:1411.4028, 2014.
  * Arute et al. [2019] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando GSL Brandao, David A Buell, et al.  Quantum supremacy using a programmable superconducting processor.  _Nature_ , 574(7779):505–510, 2019.  doi: 10.1038/s41586-019-1666-5.
  * Valiant [2002] Leslie G Valiant.  Expressiveness of matchgates.  _Theor. Comput. Sci._ , 289(1):457–471, 2002\.  doi: 10.1016/S0304-3975(01)00325-5.
  * Knill [2001] Emanuel Knill.  Fermionic linear optics and matchgates.  arXiv:quant-ph/0108033, 2001.
  * Terhal and DiVincenzo [2002] Barbara M Terhal and David P DiVincenzo.  Classical simulation of noninteracting-fermion quantum circuits.  _Phys. Rev. A_ , 65(3):032325, 2002.  doi: 10.1103/PhysRevA.65.032325.
  * DiVincenzo and Terhal [2005] D. P. DiVincenzo and B. M. Terhal.  Fermionic linear optics revisited.  _Found. Phys._ , 35(12):1967–1984, 2005.  doi: 10.1007/s10701-005-8657-0.
  * Bravyi [2005] Sergey Bravyi.  Lagrangian representation for fermionic linear optics.  _Quantum Inf. Comput._ , 5:216–238, 2005.  doi: 10.26421/qic5.3-3.
  * Jozsa and Miyake [2008] Richard Jozsa and Akimasa Miyake.  Matchgates and classical simulation of quantum circuits.  _Proc. Royal Soc. A_ , 464(2100):3089–3106, 2008\.  doi: 10.1098/rspa.2008.0189.
  * Schuch and Siewert [2003] Norbert Schuch and Jens Siewert.  Natural two-qubit gate for quantum computation using the XY interaction.  _Phys. Rev. A_ , 67(3):032301, 2003.  doi: 10.1103/PhysRevA.67.032301.
  * Abrams et al. [2020] Deanna M. Abrams, Nicolas Didier, Blake R. Johnson, Marcus P. da Silva, and Colm A. Ryan.  Implementation of XY entangling gates with a single calibrated pulse.  _Nature Electronics_ , 2020.  doi: 10.1038/s41928-020-00498-1.
  * Ishikawa and Wakayama [2006] Masao Ishikawa and Masato Wakayama.  Applications of minor summation formula III, Plücker relations, lattice paths and Pfaffian identities.  _J. Combin. Theory Ser. A_ , 113(1):113–155, 2006\.  doi: 10.1016/j.jcta.2005.05.008.
  * Diaconis and Saloff-Coste [2000] Persi Diaconis and Laurent Saloff-Coste.  Bounds for Kac’s master equation.  _Commun. Math. Phys._ , 209:729–755, 2000.  doi: 10.1007/s002200050036.
  * McKay et al. [2017] David C McKay, Christopher J Wood, Sarah Sheldon, Jerry M Chow, and Jay M Gambetta.  Efficient Z gates for quantum computing.  _Phys. Rev. A_ , 96(2):022330, 2017.  doi: 10.1103/PhysRevA.96.022330.
  * Koenig and Smolin [2014] Robert Koenig and John A Smolin.  How to efficiently select an arbitrary Clifford group element.  _J. Math. Phys._ , 55(12):122202, 2014.  doi: 10.1063/1.4903507.
  * Jiang [2017] Yunjiang Jiang.  Kac’s random walk on the special orthogonal group mixes in polynomial time.  _Proc. Amer. Math. Soc._ , 145(10):4533–4541, 2017.  doi: 10.1090/proc/13598.
  * [29] Data & code supplementary to the paper "Matchgate benchmarking: Scalable benchmarking of a continuous family of many-qubit gates".  doi: 10.5281/zenodo.5833362.
  * Jozsa et al. [2015] Richard Jozsa, Akimasa Miyake, and Sergii Strelchuk.  Jordan-Wigner formalism for arbitrary 2-input 2-output matchgates and their classical simulation.  _Quant. Inform. Comp._ , 15:0541–0556, 2015.  doi: 10.26421/qic15.7-8-1.
  * Magesan and Gambetta [2020] Easwar Magesan and Jay M Gambetta.  Effective Hamiltonian models of the cross-resonance gate.  _Phys. Rev. A_ , 101(5):052308, 2020.  doi: 10.1103/PhysRevA.101.052308.
  * Sheldon et al. [2016] Sarah Sheldon, Easwar Magesan, Jerry M Chow, and Jay M Gambetta.  Procedure for systematically tuning up cross-talk in the cross-resonance gate.  _Phys. Rev. A_ , 93(6):060302, 2016.  doi: 10.1103/PhysRevA.93.060302.
  * Brod and Childs [2014] Daniel J Brod and Andrew M Childs.  The computational power of matchgates and the XY interaction on arbitrary graphs.  _Quantum Inf. Comput._ , 14:901–916, 2014.  doi: 10.26421/qic14.11-12-1.
  * DiCarlo et al. [2009] Leonardo DiCarlo, Jerry M Chow, Jay M Gambetta, Lev S Bishop, Blake R Johnson, DI Schuster, J Majer, Alexandre Blais, Luigi Frunzio, SM Girvin, et al.  Demonstration of two-qubit algorithms with a superconducting quantum processor.  _Nature_ , 460(7252):240–244, 2009.  doi: 10.1038/nature08121.
  * Long et al. [2021] J. Long, T. Zhao, M. Bal, R. Zhao, G. S. Barron, H.-S. Ku, J. A. Howard, X. Wu, C. R. H. McRae, X.-H. Deng, et al.  A universal quantum gate set for transmon qubits with strong ZZ interactions.  arXiv:2103.12305, 2021.
  * Claes et al. [2020] J. Claes, E. Rieffel, and Z. Wang.  Character randomized benchmarking for non-multiplicity-free groups with applications to subspace, leakage, and matchgate randomized benchmarking.  2020\.  doi: 10.1103/PRXQuantum.2.010351.
  * Goodman and Wallach [2009] Roe Goodman and Nolan R Wallach.  _Symmetry, representations, and invariants_ , volume 255.  Springer, 2009.
  * Fulton and Harris [2013] William Fulton and Joe Harris.  _Representation theory: a first course_ , volume 129.  Springer, 2013.
  * Kong [2021] Linghang Kong.  A framework for randomized benchmarking over compact groups.  arXiv:2111.10357, 2021.
  * Bravyi and Gosset [2017] Sergey Bravyi and David Gosset.  Complexity of quantum impurity problems.  _Commun. Math. Phys._ , 356(2):451–500, 2017\.  doi: 10.1007/s00220-017-2976-9.

##  Appendix A Proof of the variance bound

In this section we give proofs of various technical claims made in the main
text. Consider the representation $\omega$ of the group $\operatorname{O}(2n)$
on the vector space $\Gamma$ of $n$-qubit linear operators, given by
$\omega(Q)(\rho)=U(Q)\rho U(Q)^{\dagger}$ for $Q\in\operatorname{O}(2n)$ and
$\rho\in\Gamma$. This corresponds to the conjugate action of the generalized
matchgate group. Recall that the Majorana product operators $\gamma[S]$ for
$S\subseteq[2n]$ form a basis of the space $\Gamma$, and that we defined the
subspaces $\Gamma_{k}=\langle\gamma[S]\;|\;S\subseteq[2n],\;\lvert
S\rvert=k\rangle\subseteq\Gamma$ for $k\in\\{0,1,\dots,2n\\}$. It is clear
from Eq. 1 that the subspaces $\Gamma_{k}\subseteq\Gamma$ are invariant, i.e.,
that $\omega(Q)(\Gamma_{k})\subseteq\Gamma_{k}$. Hence we can consider the
restrictions $\omega_{k}$ of $\omega$ to $\Gamma_{k}$.

######  Lemma 3.

The representation $\omega$ of $\operatorname{O}(2n)$ decomposes as a direct
sum of $2n$ irreducible subrepresentations $\omega_{k}$ for $k\in\\{0,\ldots
2n\\}$, which are all inequivalent.

###### Proof.

We only need to prove that the representations $\omega_{k}$ are irreducible
and inequivalent. Consider the linear map
$\Phi_{k}\colon\Gamma_{k}\mapsto\wedge^{k}\mathbb{C}^{2n}$ that sends each
$\gamma[S]$ for $S\subseteq[2n]$ to the antisymmetric tensor product
$\wedge_{s\in S}\lvert s\rangle$. This is an isomorphism and moreover
$\Phi_{k}(\omega(Q)\rho)=Q^{\otimes k}\Phi_{k}(\rho)$, as follows from Eq. 1.
Thus we can infer the irreducibility and mutual inequivalence of $\omega_{k}$
from the representation theory of $\operatorname{O}(2n)$ on anti-symmetric
tensor powers, which is well-known (see, e.g., [37, Cor. 5.5.6.]). ∎

We now give proof of our variance bound in the noise-free limit.

######  Theorem 1 (restated).

Consider the estimator $\hat{f}_{k}(m)$ for the quantity $f_{k}(m)$ defined in
Eqs. 2 and 5. Assuming no noise, its variance is bounded (uniformly in $k$ and
$m$) as

$\mathds{V}(\hat{f}_{k}(m))=\frac{1}{K}O(\mathrm{poly}(n)).$

###### Proof.

Begin by considering general correlation functions
$\alpha_{k}(x,Q)=N_{k}^{-1}\operatorname{Tr}(E_{x}P_{k}\omega(Q)(\rho_{0}))$.
Consider the estimator $\hat{f}_{k}(m)$ obtained by performing Algorithm 1 for
a fixed sequence length $m$, sampling $K$ random sequences and performing $L$
measurements per sequence. For any fixed $Q\in\operatorname{O}(2n)$, let
$X_{k}(Q)$ be a random variable taking value $\alpha_{k}(x,Q)$ with
probability $p(x|Q,m)$. Also let
$X_{k}^{\\{L\\}}(Q)=\frac{1}{L}\sum_{i=1}^{L}X^{i}_{k}(Q)$ denote the random
variable defined by averaging $L$ i.i.d. copies $X^{i}_{k}(Q)$ of $X_{k}(Q)$.
Finally let $Y_{k}$ be the random variable defined by drawing $Q$ uniformly at
random and taking the corresponding random variable $X_{k}^{\\{L\\}}(Q)$ (that
is, $Y_{k}=X_{k}^{\\{L\\}}(Q)$ where $Q$ is uniformly random). It is clear
that the mean of $Y_{k}$ is $f_{k}(m)$ ,and moreover the variance of the
estimator $\hat{f}_{k}$ is equal to $\frac{1}{K}\mathds{V}(Y_{k})$. The
variance of $Y_{k}$ is (by the law of total variation):

$\displaystyle\mathds{V}(Y_{k})$
$\displaystyle=\mathds{V}_{Q}\bigl{[}\mathds{E}\bigl{(}X_{k}^{\\{L\\}}(Q)\bigr{)}\bigr{]}+\mathds{E}_{Q}\bigl{[}\mathds{V}\bigl{(}X_{k}^{\\{L\\}}(Q)\bigr{)}\bigr{]},$
(7)

where the inner mean and variance are computed for arbitrary but fixed $Q$,
while the outer ones are computed with respect to the uniformly random choice
of $Q$. Now using the definitions of variance and expectation we get

$\displaystyle\mathds{V}(Y_{k})$
$\displaystyle=\int_{\operatorname{O}(2n)}dQ\bigg{(}\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q)p(x|Q,m)\bigg{)}^{2}-f_{k}(m)^{2}$
(8) $\displaystyle\hskip
30.00005pt+\frac{1}{L}\bigg{[}\int_{\operatorname{O}(2n)}dQ\Big{(}\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q)^{2}p(x|Q,m)\Big{)}-\int_{\operatorname{O}(2n)}dQ\,\Big{(}\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q)p(x|Q,m)\Big{)}^{2}\bigg{]},$
(9)

Throwing away the negative terms we get

$\mathds{V}(Y_{k})\leq\frac{L-1}{L}\int_{\operatorname{O}(2n)}dQ\,\sum_{x,x^{\prime}\in\\{0,1\\}^{n}}\alpha_{k}(x,Q)\alpha_{k}(x^{\prime},Q)p(x|Q,m)p(x^{\prime}|Q,m)+\frac{1}{L}\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q)^{2}p(x|Q,m).$
(10)

Now using the fact that the $\alpha_{k}(x,Q)$ are real functions and that for
all $a,b\in\mathbb{R}$ we have $2ab\leq a^{2}+b^{2}$ we can simplify this
further to

$\mathds{V}(Y_{k})\leq\frac{L-1}{L}2^{n}\int_{\operatorname{O}(2n)}dQ\,\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q)^{2}p(x|Q,m)^{2}+\frac{1}{L}\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q)^{2}p(x|Q,m).$
(11)

Now define $\lvert\theta_{x}^{k}\rangle=\lvert x\rangle$ for even $k$ and
$\lvert\theta_{x}^{k}\rangle=H^{\otimes n}\lvert x\rangle$ for odd $k$ (where
$H$ is the single qubit Hadamard operator).

We begin by noting that for both $k$ even and $k$ odd there always exists a
generalized matchgate $U(Q_{x}^{k})$ s.t.
$U(Q_{x}^{k})\lvert\theta_{0}^{k}\rangle=\lvert\theta_{x}^{k}\rangle$. Hence
by the invariance of the Haar measure under left multiplication and the fact
that $P_{k}$ commutes with $\omega(Q_{x}^{k})$ we have that

$\mathds{V}(Y_{k})\leq\frac{L-1}{L}2^{2n}\int_{\operatorname{O}(2n)}dQ\alpha_{k}(0,Q)^{2}p(0|Q,m)^{2}+2^{n}\frac{1}{L}\alpha_{k}(0,Q)^{2}p(0|Q,m).$
(12)

We can drop the constant factors of $L$ and upper bound the mixed integrals by
monomial integrals (by using that
$\alpha_{k}(0,Q)^{2}p(0|Q,m)^{2}\leq(\alpha_{k}(0,Q)^{4}+p(0|Q,m)^{4})/2$ :

$\displaystyle\mathds{V}(Y_{k})$ $\displaystyle\leq
2^{4n}\max\left[\int_{\operatorname{O}(2n)}dQ\,(2^{-n}\alpha_{k}(0,Q))^{4},\;\;\int_{\operatorname{O}(2n)}dQ\,p(0|Q,m)^{4}\right]$
(13) $\displaystyle\hskip
150.00023pt+2^{3n}\max\left[\int_{\operatorname{O}(2n)}dQ\,(2^{-n}\alpha_{k}(0,Q))^{3},\;\;\int_{\operatorname{O}(2n)}dQ\,p(0|Q,m)^{3}\right].$
(14)

First we consider the integrals over the correlation function $\alpha_{k}$.
Setting $t\in\\{3,4\\}$ we calculate the integral

$\displaystyle\int_{\operatorname{O}(2n)}dQ\,(2^{-n}\alpha_{k}(0,Q))^{t}$
$\displaystyle=\frac{2^{-tn}}{N_{k}^{t}}\int_{\operatorname{O}(2n)}\operatorname{Tr}\big{(}\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
P_{k}(U(Q)\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q)^{\dagger})\big{)}^{t}$ (15)
$\displaystyle=\frac{\binom{2n}{k}^{t}}{\binom{n}{\lfloor{k/2\rfloor}}^{2t}}\int_{\operatorname{O}(2n)}dQ\,\int_{\operatorname{O}(2n)}dQ^{\prime}\,\operatorname{Tr}\big{(}\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
P_{k}(U(Q^{\prime})U(Q)\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q)^{\dagger}U(Q^{\prime})^{\dagger})\big{)}^{t}$ (16)
$\displaystyle=\frac{\binom{2n}{k}^{t}}{\binom{n}{\lfloor{k/2\rfloor}}^{2t}}\operatorname{Tr}\bigg{(}\int_{\operatorname{O}(2n)}dQ^{\prime}\,U(Q^{\prime})\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q^{\prime})^{\dagger}P_{k}\bigg{(}\int_{\operatorname{O}(2n)}dQ\,U(Q)\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q)^{\dagger}\bigg{)}\bigg{)}^{t},$ (17)

using the invariance of the Haar measure and the fact that $P_{k}$ commutes
with the conjugate action $\omega(Q)$. Now, since $P_{k}$ is an orthogonal
projector, it is a contraction in the Hilbert-Schmidt norm, and we get (using
Cauchy-Schwartz)

$\displaystyle\int_{\operatorname{O}(2n)}dQ\,(2^{-n}\alpha_{k}(0,Q))^{t}$
$\displaystyle\leq\frac{\binom{2n}{k}^{t}}{\binom{n}{\lfloor{k/2\rfloor}}^{2t}}\operatorname{Tr}\bigg{(}\int_{\operatorname{O}(2n)}dQ^{\prime}U(Q^{\prime})\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q^{\prime})^{\dagger}\int_{\operatorname{O}(2n)}dQ\,U(Q)\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q)^{\dagger}\bigg{)}^{t}.$ (18)

Using the invariance of the Haar measure (to absorb one of the integrals) and
the definition of $p(0|Q,m)$ we see that the RHS becomes

$\displaystyle\frac{\binom{2n}{k}^{t}}{\binom{n}{\lfloor{k/2\rfloor}}^{2t}}$
$\displaystyle\operatorname{Tr}\bigg{(}\int_{\operatorname{O}(2n)}dQ^{\prime}U(Q^{\prime})\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q^{\prime})^{\dagger}\int_{\operatorname{O}(2n)}dQ\,U(Q)\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q)^{\dagger}\bigg{)}^{t}$ (19)
$\displaystyle=\frac{\binom{2n}{k}^{t}}{\binom{n}{\lfloor{k/2\rfloor}}^{2t}}\int_{\operatorname{O}(2n)}dQ\,\operatorname{Tr}\big{(}\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert)U(Q)\lvert\theta^{k}_{0}\rangle\\!\langle\theta^{k}_{0}\rvert
U(Q)^{\dagger}\big{)}^{t}$ (20)
$\displaystyle=\frac{\binom{2n}{k}^{t}}{\binom{n}{\lfloor{k/2\rfloor}}^{2t}}\int_{\operatorname{O}(2n)}dQ\,p(0|Q,m)^{t},$
(21)

and hence

$\displaystyle\int_{\operatorname{O}(2n)}dQ\,(2^{-n}\alpha_{k}(0,Q))^{t}\leq\frac{\binom{2n}{k}^{t}}{\binom{n}{\lfloor{k/2\rfloor}}^{2t}}\int_{\operatorname{O}(2n)}dQ\,p(0|Q,m)^{t}.$
(22)

Now we use Stirling’s approximation ($\sqrt{2\pi
n}(n/e)^{n}e^{\frac{1}{12n+1}}\leq n!\leq\sqrt{2\pi
n}(n/e)^{n}e^{\frac{1}{12n}}$) to note that
$\frac{\binom{2n}{k}}{\binom{n}{\lfloor{k/2\rfloor}}^{2}}=O\big{(}\mathrm{poly}(n,k)\big{)}$
we thus only need to consider the integral over $p(0|Q,m)$:

$\displaystyle\int_{\operatorname{O}(2n)}dQ\,p(0|Q,m)^{t}$
$\displaystyle=\int_{\operatorname{O}(2n)}dQ\operatorname{Tr}\big{(}\lvert\theta_{0}^{k}\rangle\\!\langle\theta_{0}^{k}\rvert
U(Q)\lvert\theta_{0}^{k}\rangle\\!\langle\theta_{0}^{k}\rvert
U(Q)^{\dagger}\big{)}^{t}$ (23)
$\displaystyle=\int_{\operatorname{SO}(2n)}dR\operatorname{Tr}\big{(}\lvert\theta_{0}^{k}\rangle\\!\langle\theta_{0}^{k}\rvert
U(R)\lvert\theta_{0}^{k}\rangle\\!\langle\theta_{0}^{k}\rvert
U(R)^{\dagger}\big{)}^{t}$ (24) $\displaystyle\hskip
100.00015pt+\int_{\operatorname{SO}(2n)}dR\operatorname{Tr}\big{(}\lvert\theta_{0}^{k}\rangle\\!\langle\theta_{0}^{k}\rvert
U(R)X_{n}\lvert\theta_{0}^{k}\rangle\\!\langle\theta_{0}^{k}\rvert
X_{n}U(R)^{\dagger}\big{)}^{t},$ (25)

where $X_{n}$ is a bit flip on the last qubit. Defining $\lvert
e\rangle=X_{n}\lvert 0\rangle^{\otimes n}$ and noting that
$X_{n}\lvert+\rangle^{\otimes n}=\lvert+\rangle^{\otimes n}$ we see that for
both $k$ even and $k$ odd:

$\int_{\operatorname{O}(2n)}dQ\,p(0|Q,m)^{t}\leq
2\max\bigg{\\{}\int_{\operatorname{SO}(2n)}dR\operatorname{Tr}\big{(}\lvert\phi\rangle\\!\langle\phi\rvert
U(R)\lvert\phi\rangle\\!\langle\phi\rvert
U(R)^{\dagger}\big{)}^{t}\;\;\bigg{|}\;\;\lvert\phi\rangle\in\\{\lvert
0\rangle,\lvert e\rangle,\frac{1}{\sqrt{2}}(\lvert 0\rangle+\lvert
e\rangle)\\}\bigg{\\}},$ (26)

where we used that $\lvert+\rangle^{\otimes
n}=U(R_{+})\frac{1}{\sqrt{2}}(\lvert 0\rangle+\lvert e\rangle)$ for some
$R_{+}\in\operatorname{SO}(2n)$ as well as Haar invariance. All these
integrals are bounded by the moment bounds in Theorem 6, and we find that

$\int_{\operatorname{O}(2n)}dQ\,p(0|Q,m)^{t}=2^{-tn}\,O(\mathrm{poly}(n)).$
(27)

Hence we have

$\displaystyle\mathds{V}(Y_{k})=O(\mathrm{poly}(n)),$ (28)

for all $k\in[2n]$, which proves the theorem. ∎

##  Appendix B Proof of the moment bound

In this section we state and prove Theorem 6, which contains our moment bounds
and implies Lemma 2 in the main text. The $n$-qubit Hilbert space
$(\mathbb{C}^{2})^{\otimes n}$ carries a representation of the Lie algebra
$\mathfrak{so}(2n)$. The corresponding Lie group representation corresponds
precisely to the matchgate action, so that Eq. 1 holds. We now recall some
notions from representation theory (see, e.g., [38] for a gentle
introduction). Define fermionic creation and annihilation operators by

$\displaystyle a_{i}=\frac{\gamma_{2i-1}+i\gamma_{2i}}{2},\quad
a_{i}^{\dagger}=\frac{\gamma_{2i-1}-i\gamma_{2i}}{2}\quad\text{ for }i\in[n],$
(29)

in terms of the Majorana fermion operators, which were defined in the main
text as operators on $(\mathbb{C}^{2})^{\otimes n}$. In this way, we can
identify $(\mathbb{C}^{2})^{\otimes n}$ with the fermionic Fock space
$\bigwedge\mathbb{C}^{n}$. Note that all-zero basis state $\lvert 0\rangle$
corresponds to the fermionic Fock vacuum, since we have $a_{i}\lvert
0\rangle=0$ for $i\in[n]$. Now, consider the following operators:

$\displaystyle
X^{\prime}_{ij}=a_{i}^{\dagger}a_{j}-\frac{1}{2}\delta_{ij},\quad
Y^{\prime}_{ij}=a_{i}^{\dagger}a_{j}^{\dagger}\;\;\text{(for $i<j$)},\quad
Z^{\prime}_{ij}=a_{i}a_{j}\;\;\text{(for $i<j$)}.$ (30)

Here we follow the notation and conventions of [38] (and caution that these
operators are not the Pauli matrices). The operators defined in Eq. 30 satisfy
the commutation relations of $\mathfrak{so}(2n)$, and hence define a
representation of $\mathfrak{so}(2n)$ on the $n$-qubit Hilbert space. This
representation generates the matchgate group.

Any irreducible representation of $\mathfrak{so}(2n)$ is classified by its
highest weight. A _weight vector_ is a joint eigenvector of the operators
$H^{\prime}_{i}:=X^{\prime}_{ii}$ for $i\in[n]$; and the vector of eigenvalues
is simply called the _weight_. Note that, as an operator on the $n$-qubit
Hilbert space, $H^{\prime}_{i}$ is nothing but $-\frac{1}{2}Z_{i}$ ($Z_{i}$ is
the Pauli $Z$-matrix acting on the $i$-th qubit). Moreover, a _highest weight
vector_ is a weight vector that is annihilated by the operators
$\\{X^{\prime}_{ij}\\}_{i<j}$ and $\\{Y^{\prime}_{ij}\\}_{i<j}$. Any
irreducible representation contains a unique (up to phase) highest weight
vector, and its weight, called the _highest weight_ , characterizes the
representation completely. It is well-known that the representation of
$\mathfrak{so}(2n)$ on $(\mathbb{C}^{2})^{\otimes n}$ defined above decomposes
into two irreducible representations, with highest weights

$\displaystyle\alpha:=(\tfrac{1}{2},\dots,\tfrac{1}{2},\tfrac{1}{2})\quad\text{and}\quad\beta:=(\tfrac{1}{2},\dots,\tfrac{1}{2},-\tfrac{1}{2})\in(\mathbb{Z}/2)^{n}.$
(31)

If $n$ is even, the subrepresentation with highest weight $\alpha$ is the even
particle number subspace of $(\mathbb{C}^{2})^{\otimes n}$, and the
subrepresentation with highest weight $\beta$ is the odd particle number
subspace. If $n$ is odd, then the reverse is true: $\beta$ corresponds to the
even particle number subspace, and $\alpha$ corresponds to the odd particle
number subspace. See [38, Prop. 20.15]. Indeed, the all-one state $\lvert
1\rangle^{\otimes n}$ is a highest weight vector of weight $\alpha$ (but its
parity depends on the parity of $n$), while the state $X_{n}\lvert
1\rangle^{\otimes n}$ is a highest weight vector of weight $\beta$.

In our analysis, it will be useful to instead consider the _lowest weight
vectors_ , which are the weight vectors that are annihilated by the
$\\{X^{\prime}_{ij}\\}_{i>j}$ and $\\{Z^{\prime}_{ij}\\}_{i<j}$. Just like the
highest weight vectors, they characterize the irreducible representation
uniquely. Clearly, both the all-zero state $\lvert 0\rangle$ and the state
$\lvert e\rangle=a_{n}^{\dagger}\lvert 0\rangle=X_{n}\lvert 0\rangle$ are
lowest weight vectors. We can compute the _highest_ weight of the
corresponding representation by observing that $\lvert 0\rangle$ is the vacuum
(an even particle number state), while $\lvert e\rangle$ is a single-particle
state (an odd particle number state). We summarize:

######  Lemma 4.

The vectors $\lvert 0\rangle$ and $\lvert e\rangle=a^{\dagger}_{n}\lvert
0\rangle=X_{n}\lvert 0\rangle$ are lowest weight vectors in
$(\mathbb{C}^{2})^{\otimes n}$. For even $n$, $\lvert 0\rangle$ is contained
in an irreducible subrepresentation with highest weight $\alpha$ and $\lvert
e\rangle$ is contained in an irreducible subrepresentation with highest weight
$\beta$. For odd $n$, the opposite holds: $\lvert 0\rangle$ corresponds to
highest weight $\beta$ and $\lvert 1\rangle$ to highest weight $\alpha$.

In general, the highest weight can be obtained from the lowest weight by the
action of the Weyl group, which can permute the entries of the weight as well
as swap an even number of signs, until we obtain a weight
$\omega\in(\mathbb{Z}/2)^{n}$ that satisfies
$\omega_{1}\geq\dots\geq\omega_{n-1}\geq\lvert\omega_{n}\rvert$. For $\lvert
0\rangle$, the weight is $(-\tfrac{1}{2},\dots,-\frac{1}{2},-\frac{1}{2})$, so
we obtain $\alpha$ if $n$ is even and $\beta$ if $n$ is odd. For $\lvert
e\rangle$, the weight is $(-\tfrac{1}{2},\dots,-\frac{1}{2},\frac{1}{2})$, so
we obtain $\beta$ if $n$ is even and $\alpha$ if $n$ is odd.

It is well-known and easy to see that the tensor product of highest weight
vectors is again a highest weight vector, with associated highest weight the
sum of the highest weights of the individual tensor factors (see [38, Obs.
13.2]). The same is true for lowest weight vectors. Accordingly, for any $t$
and $m\in[t]$, we can define the following lowest weight vector:

$\displaystyle\lvert\Omega^{t}_{m}\rangle:=\lvert 0\rangle^{\otimes
m}\otimes\lvert e\rangle^{\otimes t-m}.$ (32)

The corresponding highest weight $\lambda_{m}^{t}\in(\mathbb{Z}/2)^{n}$,
computed as described above, is the following:

$\lambda_{m}^{t}:=\left\\{\begin{array}[]{cc}m\alpha+(t-m)\beta=(\frac{t}{2},\cdots,\frac{t}{2},m-\frac{t}{2})&\text{
if $n$ is even,}\\\
m\beta+(t-m)\alpha=(\frac{t}{2},\cdots,\frac{t}{2},\frac{t}{2}-m)&\text{ if
$n$ is odd.}\end{array}\right.$ (33)

These arguments will crucially feature in our moment bound. Indeed, the
computation of the $2t$-th moment can be reduced to an integral involving the
vectors $\lvert\Omega^{t}_{m}\rangle$ for $m\in[t]$. Since each
$\lvert\Omega^{t}_{m}\rangle$ is a lowest weight vector, it is supported in a
_single_ irreducible representation. Therefore, we can use powerful tools from
representation theory such as Schur’s lemma and the Weyl dimension formula to
compute the corresponding integrals.

Before we can state our moment calculations, we must discuss one subtlety. Not
every representation of the Lie algebra $\mathfrak{so}(2n)$ integrates to a
representation of the Lie group $\operatorname{SO}(2n)$. In particular, this
problem occurs for the representation on $\mathfrak{so}(2n)$ on
$(\mathbb{C}^{2})^{\otimes n}$ discussed above – meaning that the unitaries
$U(R)$ for $R\in\operatorname{SO}(2n)$ do _not_ define a representation of
$\operatorname{SO}(2n)$. So far, this was not important for our analysis,
since the conjugation with $U(R)$ as in Eq. 1 gives a well-defined
representation of $\operatorname{SO}(2n)$. However, it will be necessary to be
mindful of this subtlety in what follows. Fortunately, for any representation
of $\mathfrak{so}(2n)$ we always have a corresponding representation of the
Lie group $\operatorname{Spin}(2n)$, which has the same Lie algebra but is a
simply connected double cover of $\operatorname{SO}(2n)$. We denote the
$\operatorname{Spin}(2n)$-representation obtained in this way from the
$\mathfrak{so}(2n)$-representation on $(\mathbb{C}^{2})^{\otimes n}$ by
$\rho(R)$ for $R\in\operatorname{Spin}(2n)$. The matchgate group
$\mathcal{M}_{n}$ is nothing but the image of $\operatorname{Spin}(2n)$ under
this representation.

######  Lemma 5.

For integers $t$ and complex numbers $x,y$ we have

$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,\left|\langle
0\rvert(\overline{x}+\overline{y}a_{n})U(R)(x+ya_{n}^{\dagger})\lvert
0\rangle\right|^{2t}$ $\displaystyle\hskip
10.00002pt=(n-1)!^{2}\sum_{m=0}^{t}|x|^{4m}|y|^{4(t-m)}\binom{t}{m}^{2}\prod_{i=1}^{t}\frac{(i+n-1)!}{\sqrt{i!(2n-2+i)!}}\sqrt{\binom{n+t/2-1}{t/2}}\frac{m!}{(m+n-1)!}\frac{(t-m)!}{(t-m+n-1)!}.$
(34)

###### Proof.

We first write this as a Haar integral over the spin group, which does not
change the value as explained above.

$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,\left|\langle
0\rvert(\overline{x}+\overline{y}a_{n})U(R)(x+ya_{n}^{\dagger})\lvert
0\rangle\right|^{2t}=\int_{\operatorname{Spin}(2n)}dR\,\,\left|\langle
0\rvert(\overline{x}+\overline{y}a_{n})\rho(R)(x+ya_{n}^{\dagger})\lvert
0\rangle\right|^{2t}.$ (35)

Recall that $\lvert e\rangle=a_{n}^{\dagger}\lvert 0\rangle$. The Gaussian
fermionic unitaries preserve the parity of fermions, therefore,

$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,$ $\displaystyle\left|\langle
0\rvert(\overline{x}+\overline{y}a_{n})\rho(R)(x+ya_{n}^{\dagger})\lvert
0\rangle\right|^{2t}$ (36)
$\displaystyle=\int_{\operatorname{SO}(2n)}dR\,\,\left||x|^{2}\langle
0\rvert\rho(R)\lvert 0\rangle+|y|^{2}\langle e\rvert\rho(R)\lvert
e\rangle\right|^{2t}$ (37)
$\displaystyle=\int_{\operatorname{SO}(2n)}dR\,\,\left|\sum_{m=0}^{t}\binom{t}{m}|x|^{2m}|y|^{2(t-m)}\langle
0\rvert\rho(R)\lvert 0\rangle^{m}\langle e\rvert\rho(R)\lvert
e\rangle^{t-m}\right|^{2}$ (38)
$\displaystyle=\\!\\!\sum_{m,m^{\prime}=0}^{t,t}\binom{t}{m}\binom{t}{m^{\prime}}|x|^{2m+2m^{\prime}}|y|^{2(t-m)+2(t-m^{\prime})}\bigg{(}\int_{\operatorname{SO}(2n)}dR\,\langle
0\rvert\rho(R)\lvert 0\rangle^{m}\langle 0\rvert\rho(R)^{\dagger}\lvert
0\rangle^{m^{\prime}}$ (39) $\displaystyle\hskip 270.00041pt\times\langle
e\rvert\rho(R)\lvert e\rangle^{t-m}\langle e\rvert\rho(R)^{\dagger}\lvert
e\rangle^{t-m^{\prime}}\bigg{)},$ (40)

where the first equality follows from preserving the parity. The rest are
simple algebraic manipulations. Now recall that we defined
$\lvert\Omega^{t}_{m}\rangle:=\lvert 0\rangle^{\otimes m}\otimes\lvert
e\rangle^{\otimes t-m}$ in Eq. 32. Thus, Eq. 36 can be written as

$\displaystyle\int_{\operatorname{SO}(2n)}$ $\displaystyle dR\,\,\left|\langle
0\rvert(\overline{x}+\overline{y}a_{n})\rho(R)(x+ya_{n}^{\dagger})\lvert
0\rangle\right|^{2t}$
$\displaystyle=\\!\\!\\!\\!\sum_{m,m^{\prime}=0}^{t,t}\binom{t}{m}\binom{t}{m^{\prime}}|x|^{2m+2m^{\prime}}|y|^{2(t-m)+2(t-m^{\prime})}\left(\int_{\operatorname{Spin}(2n)}dR\,\langle\Omega^{t}_{m}\rvert\rho^{\otimes
t}(R)\lvert\Omega^{t}_{m}\rangle\langle\Omega^{t}_{m^{\prime}}\rvert\rho^{\otimes
t}(R)^{\dagger}\lvert\Omega^{t}_{m^{\prime}}\rangle\right).$ (41)

Now, because $\lvert\Omega^{t}_{m}\rangle$ is a lowest weight vector, its
support is limited to the irreducible representation of the $\text{so}(2n)$
Lie algebra given by the corresponding highest weight $\lambda_{m}^{t}$ as
given in Eq. 33. We denote this irreducible representation by
$\rho_{\lambda_{m}^{t}}$. Therefore, the above relation can be re-written as

$\displaystyle\int_{\operatorname{SO}(2n)}$ $\displaystyle dR\,\,\left|\langle
0\rvert(\overline{x}+\overline{y}a_{n})\rho(R)(x+ya_{n}^{\dagger})\lvert
0\rangle\right|^{2t}$
$\displaystyle=\\!\\!\\!\\!\sum_{m,m^{\prime}=0}^{t,t}\binom{t}{m}\binom{t}{m^{\prime}}|x|^{2m+2m^{\prime}}|y|^{2(t-m)+2(t-m^{\prime})}\left(\int_{\operatorname{Spin}(2n)}dR\,\langle\Omega^{t}_{m}\rvert\rho_{\lambda_{m}^{t}}(R)\lvert\Omega^{t}_{m}\rangle\langle\Omega^{t}_{m^{\prime}}\rvert\rho_{\lambda_{m^{\prime}}}(R)^{\dagger}\lvert\Omega^{t}_{m^{\prime}}\rangle\right).$
(42)

Let us focus on the integral in the parenthesis. As a result of the Schur
orthogonality relations, it is straightforward to see that the Haar integral
$\int{dR\rho_{\lambda_{m}^{t}}(R)\lvert\Omega^{t}_{m}\rangle\langle\Omega^{t}_{m^{\prime}}\rvert\rho_{\lambda_{m^{\prime}}}(R)^{\dagger}}$
vanishes if $m\neq m^{\prime}$. Furthermore, if $m=m^{\prime}$, then

$\int_{\operatorname{Spin}(2n)}{dR\,\rho_{\lambda_{m}^{t}}(R)\lvert\Omega^{t}_{m}\rangle\langle\Omega^{t}_{m}\rvert\rho_{\lambda_{m}}(R)^{\dagger}}\propto\rho_{\lambda_{m}}(\text{Id}),$
(43)

as a consequence of Schur’s lemma. The proportionality constant can be
calculated by comparing the traces of the both side of the equality, and it is
equal to $1/\dim({\rho_{\lambda_{m}^{t}}})$. Hence, we obtain,

$\int_{\operatorname{Spin}(2n)}dR\,\,\left|\langle
0\rvert(\overline{x}+\overline{y}a_{n})\rho(R)(x+ya_{n}^{\dagger})\lvert
0\rangle\right|^{2t}=\sum_{m=0}^{t}\binom{t}{m}^{2}|x|^{4m}|y|^{4(t-m)}\frac{1}{\dim({\rho_{\lambda_{m}^{t}}})}.$
(44)

The dimension of $\rho_{\lambda_{m}^{t}}$ can be directly calculated using
explicit relations for the dimension of the irreducible representations of
$\mathfrak{so}(2n)$. For an arbitrary highest weight
$\mu=(\mu_{1},\mu_{2},\cdots,\mu_{n})$, the dimension of the corresponding
irreducible representation is given by one of Weyl dimension formulas [38, Eq.
24.41]:

$\dim{(\rho_{\mu})}=\prod_{1\leq i<j\leq
n}\frac{l_{i}^{2}-l_{j}^{2}}{m_{i}^{2}-m_{j}^{2}},\quad\text{with
}m_{i}:=n-i,\text{ and }l_{i}:=\mu_{i}+n-i.$

After a few lines of algebra, this leads to,

$\dim{(\rho_{\lambda_{m}^{t}})}=\left(\prod_{1\leq i<j\leq
n-1}\frac{t+i+j}{i+j}\right)\times\frac{1}{(n-1)!^{2}}\times\frac{(m+n-1)!}{m!}\times\frac{(t-m+n-1)!}{(t-m)!}.$
(45)

We focus on the term inside the parentheses. By some manipulation of the
product factors and the definition of the binomial we see:

$\displaystyle\prod_{1\leq i<j\leq n-1}\frac{t+i+j}{i+j}$
$\displaystyle=\sqrt{\prod_{1\leq i<j\leq
n-1}\frac{t+i+j}{i+j}\times\prod_{1\leq j<i\leq n-1}\frac{t+i+j}{i+j}}$ (46)
$\displaystyle=\sqrt{\prod_{\begin{subarray}{c}i,j=1\\\ i\neq
j\end{subarray}}^{n-1}\frac{t+i+j}{i+j}}$ (47)
$\displaystyle=\sqrt{\prod_{i=1}^{n-1}\frac{i}{i+\frac{t}{2}}}\times\sqrt{\prod_{i,j=1}^{n-1}\frac{t+i+j}{i+j}}$
(48)
$\displaystyle=\binom{\frac{t}{2}+n-1}{\frac{t}{2}}^{-1/2}\times\sqrt{\prod_{i,j=1}^{n-1}\frac{t+i+j}{i+j}}.$
(49)

Furthermore, we see that

$\displaystyle\prod_{i,j=1}^{n-1}\frac{t+i+j}{i+j}$
$\displaystyle=\prod_{i=1}^{n-1}\left(\frac{\prod_{j=1}^{n-t-1}(t+i+j)}{\prod_{j=t+1}^{n-1}(i+j)}\times\frac{\prod_{j=n-t}^{n-1}(t+i+j)}{\prod_{j=1}^{t}(i+j)}\right)$
(50)

which is nothing more than splitting the products over the index $j$ into two
components set by $t$. Working this out further we get

$\displaystyle\prod_{i,j=1}^{n-1}\frac{t+i+j}{i+j}$
$\displaystyle=\prod_{i=1}^{n-1}\left(\frac{\prod_{j=1}^{n-t-1}(t+i+j)}{\prod_{j=1}^{n-1-t}(t+i+j)}\times\frac{\prod_{j=1}^{t}(n-1+i+j)}{\prod_{j=1}^{t}(i+j)}\right)$
(51)
$\displaystyle=\prod_{i=1}^{n-1}\left(\frac{\prod_{j=1}^{t}(n-1+i+j)}{\prod_{j=1}^{t}(i+j)}\right)$
(52)
$\displaystyle=\prod_{i=1}^{t}\left(\frac{\prod_{j=1}^{n-1}(n-1+i+j)}{\prod_{j=1}^{n-1}(i+j)}\right)$
(53) $\displaystyle=\prod_{i=1}^{t}i!\frac{(2n-2+i)!}{(i+n-1)!^{2}}.$ (54)

Combining all of these relations, we have:

$\displaystyle\int_{\operatorname{Spin}(2n)}$ $\displaystyle
dR\,\,\left|\langle
0\rvert(\overline{x}+\overline{y}a_{n})\rho(R)(x+ya_{n}^{\dagger})\lvert
0\rangle\right|^{2t}$ (55)
$\displaystyle=(n-1)!^{2}\sum_{m=0}^{t}|x|^{4m}|y|^{4(t-m)}\binom{t}{m}^{2}\prod_{i=1}^{t}\frac{(i+n-1)!}{\sqrt{i!(2n-2+i)!}}\sqrt{\binom{n+t/2-1}{t/2}}\frac{m!}{(m+n-1)!}\frac{(t-m)!}{(t-m+n-1)!}.$

∎

Lastly, we prove the main result of this section.

######  Theorem 6.

Let $\lvert 0\rangle$ be the all-zero state and define the state $\lvert
e\rangle=X_{n}\lvert 0\rangle=a_{n}^{\dagger}\lvert 0\rangle$. We have, for
any fixed $t$:

$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,\left|\frac{\langle
0\rvert+\langle e\rvert}{\sqrt{2}}U(R)\frac{\lvert 0\rangle+\lvert
e\rangle}{\sqrt{2}}\right|^{2t}=2^{-tn}\,O(\mathrm{poly}(n)),$ (56)
$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,\left|\langle 0\rvert
U(R)\lvert 0\rangle\right|^{2t}=2^{-tn}\,O(\mathrm{poly}(n)),$ (57)
$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,\left|\langle e\rvert
U(R)\lvert e\rangle\right|^{2t}=2^{-tn}\,O(\mathrm{poly}(n)).$ (58)

###### Proof.

We note that the first integral is of the form Lemma 5 with
$x=\frac{1}{\sqrt{2}}$ and $y=\frac{1}{\sqrt{2}}$:

$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,$
$\displaystyle\left|\frac{\langle 0\rvert+\langle
e\rvert}{\sqrt{2}}U(R)\frac{\lvert 0\rangle+\lvert
e\rangle}{\sqrt{2}}\right|^{2t}$ (59)
$\displaystyle=4^{-t}(n-1)!^{2}\sum_{m=0}^{t}\binom{t}{m}^{2}\prod_{i=1}^{t}\frac{(i+n-1)!}{\sqrt{i!(2n-2+i)!}}\sqrt{\binom{n+t/2-1}{t/2}}\frac{m!}{(m+n-1)!}\frac{(t-m)!}{(t-m+n-1)!}$
(60)
$\displaystyle=4^{-t}(n-1)!^{2}t!^{2}\prod_{i=1}^{t}\frac{(i+n-1)!}{\sqrt{i!(2n-2+i)!}}\sqrt{\binom{n+t/2-1}{t/2}}\sum_{m=0}^{t}\frac{1}{m!(t-m)!}\frac{1}{(m+n-1)!}\frac{1}{(t-m+n-1)!}$
$\displaystyle=\frac{4^{-t}(n-1)!^{2}t!^{2}}{(t+n-1)!^{2}}\prod_{i=1}^{t}\frac{(i+n-1)!}{\sqrt{i!(2n-2+i)!}}\sqrt{\binom{n+t/2-1}{t/2}}\sum_{m=0}^{t}\binom{t+n-1}{m}\binom{t+n-1}{t-m}.$
(61)

Now, we can use the binomial identity

$\sum_{m=0}^{t}\binom{t+n-1}{m}\binom{t+n-1}{t-m}=\binom{2(t+n-1)}{t}.$ (62)

Hence, we have

$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,$
$\displaystyle\left|\frac{\langle 0\rvert+\langle
e\rvert}{\sqrt{2}}U(R)\frac{\lvert 0\rangle+\lvert
e\rangle}{\sqrt{2}}\right|^{2t}$ (63)
$\displaystyle=4^{-t}\prod_{i=1}^{t}\frac{(i+n-1)!}{\sqrt{i!(2n-2+i)!}}\sqrt{\binom{n+t/2-1}{t/2}}\binom{2(t+n-1)}{t}\binom{t+n-1}{t}^{-2}$
(64)
$\displaystyle=4^{-t}\left(\prod_{i=1}^{t}\binom{2n+2i-2}{n+i-1}^{-1/2}\right)\left(\prod_{i=1}^{t}\binom{2n+2i-2}{i}^{1/2}\right)\sqrt{\binom{n+t/2-1}{t/2}}\binom{2(t+n-1)}{t}\binom{t+n-1}{t}^{-2}.$
(65)

The last four terms in Eq. 65 (up to square roots) are polynomials in $n$.
Therefore, we have

$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,\left|\frac{\langle
0\rvert+\langle e\rvert}{\sqrt{2}}U(R)\frac{\lvert 0\rangle+\lvert
e\rangle}{\sqrt{2}}\right|^{2t}=4^{-t}\left(\prod_{i=1}^{t}\binom{2n+2i-2}{n+i-1}^{-1/2}\right)\times
O(\mathrm{poly}(n)).$ (66)

Using Stirling’s approximation $\sqrt{2\pi n}(n/e)^{n}e^{\frac{1}{12n+1}}\leq
n!\leq\sqrt{2\pi n}(n/e)^{n}e^{\frac{1}{12n}}$, the following bound holds,

$\frac{4^{n}}{\sqrt{\pi
n}}\exp\left[-\frac{1}{8n}\right]\leq\frac{4^{n}}{\sqrt{\pi
n}}\exp\left[\frac{1}{24n+1}\\!-\\!\frac{1}{6n}\right]\leq\binom{2n}{n}.$ (67)

We can apply this lower bound to every factor of $\binom{2n+2i-2}{n+i-1}$ in
the product factor of Eq. 66 to obtain an upper bound of this factor:

$\displaystyle 4^{-t}$
$\displaystyle\left(\prod_{i=1}^{t}\binom{2n+2i-2}{n+i-1}^{-1/2}\right)\leq
2^{-tn}2^{-t^{2}-2t}\pi^{t/4}\left(\frac{(n+t-1)!}{(n-1)!}\right)^{1/4}\exp\left(\frac{t}{16n}\right).$
(68)

Inserting this relation back into Eq. 66, we conclude that

$\int_{\operatorname{SO}(2n)}dR\,\,\left|\frac{\langle 0\rvert+\langle
e\rvert}{\sqrt{2}}U(R)\frac{\lvert 0\rangle+\lvert
e\rangle}{\sqrt{2}}\right|^{2t}=2^{-tn}\times O(\mathrm{poly}(n)).$ (69)

Next, we discuss Eqs. 57 and 58. Using Lemma 5 for $x=0,y=1$ and $x=1,y=0$ we
have

$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,\left|\langle 0\rvert
U(R)\lvert 0\rangle\right|^{2t}=$
$\displaystyle\int_{\operatorname{SO}(2n)}dR\,\,\left|\langle e\rvert
U(R)\lvert e\rangle\right|^{2t}$ $\displaystyle=$
$\displaystyle\left(\prod_{i=1}^{t}\frac{(i+n-1)!}{\sqrt{i!(2n-2+i)!}}\right)\sqrt{\binom{n+t/2-1}{t/2}}\frac{t!}{(t+n-1)!}$
$\displaystyle=$
$\displaystyle\left(\binom{2n-2}{n-1}^{-1/2}\right)^{t}\times\left(\prod_{i=1}^{t}\binom{i+n-1}{i}\binom{2n-2+i}{2n-2}^{-1/2}\right)$
$\displaystyle=$ $\displaystyle\left(\binom{2n-2}{n-1}^{-1/2}\right)^{t}\times
O(\mathrm{poly}(n)).$ (70)

Incorporating Eq. 67 into Appendix B we immediately obtain the desired results
Eqs. 57 and 58. ∎

We note that Lemma 2 follows from Theorem 6. For the all-plus state, this
follows by Haar invariance and the fact that $\lvert+\rangle^{\otimes
n}=U(R_{+})\frac{1}{\sqrt{2}}(\lvert 0\rangle+\lvert e\rangle)$ for some
$R_{+}\in\operatorname{SO}(2n)$.

##  Appendix C Derivation of decay model in the gate-independent noise case

In this section we provide an “artisanal” derivation of Eq. 4 under the
assumption of gate-independent noise. We assume that there exists a quantum
channel $\Lambda$ such that every generalized matchgate $U(Q)$ is implemented
on the device as $\Lambda\circ\omega(Q)$, where
$\omega(Q)(\rho)=U(Q)\rho\,U(Q)^{\dagger}$ as defined above. We note that this
assumption is not very realistic and the decay model can be derived under much
weaker conditions using the general arguments given in [10] (specifically
Theorem $9$ therein). However this derivation has the benefit of making it
more clear what all the moving parts are. We will give an argument that the
gate-dependent noise assumption can be relaxed in the next section.

Throughout the derivation we will make use of the matrix-transfer
representation, writing matrices $\rho,E$ as vectors
$\lvert\rho\rangle,\langle E\rvert$ with trace inner product $\langle
E|\rho\rangle=\operatorname{Tr}(E^{\dagger}\rho)$. Correspondingly,
superoperators $\Lambda$ get mapped to matrices acting as
$\Lambda\lvert\rho\rangle=\lvert\Lambda(\rho)\rangle$. This representation
maps composition to matrix multiplication and interplays correctly with tensor
products. In this picture we can write the output of the matchgate
benchmarking protocol as

$\displaystyle f_{k}(m)$
$\displaystyle=\int_{\operatorname{O}(2n)}dQ_{1}\cdots
Q_{m}\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q_{1}\cdots
Q_{m})\langle\tilde{E}_{x}\rvert\Lambda\omega(Q_{m})\cdots\Lambda\omega(Q_{1})\lvert\tilde{\rho}_{0}\rangle$
(71) $\displaystyle=N_{k}^{-1}\int_{\operatorname{O}(2n)}dQ_{1}\cdots
Q_{m}\sum_{x\in\\{0,1\\}^{n}}\langle
E_{x}\otimes\tilde{E}_{x}\rvert(P_{k}\otimes\Lambda)\omega(Q_{m})^{\otimes
2}(\mathbb{1}\otimes\Lambda)\cdots(\mathbb{1}\otimes\Lambda)\omega(Q_{1})^{\otimes
2}\lvert\rho_{0}\otimes\tilde{\rho}_{0}\rangle,$ (72)

with $\tilde{E}_{x},\tilde{\rho}_{0}$ noisy versions of their ideal
counterparts. Using standard representation theory and the decomposition of
Lemma 3 we can express the integral
$\int_{\operatorname{O}(2n)}dQ\omega(Q)^{\otimes 2}$ as a projector onto a
space spanned by the vectors

$\\{\lvert\mathrm{v}(P_{k^{\prime}})\rangle\;\;\;|\;\;k\in\\{0,\ldots,2n\\}\\},$
(73)

where

$\lvert\mathrm{v}(P_{k^{\prime}})\rangle:=2^{-n}\sum_{S\subset[2n],|S|=k^{\prime}}\lvert\gamma[S]\otimes\gamma[S]\rangle,$
(74)

with
$P_{k^{\prime}}=2^{-n}\sum_{S\subset[2n],|S|=k^{\prime}}\lvert\gamma[S]\rangle\langle\gamma[S]\rvert$
the projector onto the subrepresentation $\Gamma_{k}$. Using linearity, and
writing $|P_{k^{\prime}}|$ for the dimension of $\Gamma_{k^{\prime}}$, we can
insert this into the expression for $f_{k}(m)$

$\displaystyle f_{k}(m)$
$\displaystyle=N_{k}^{-1}\sum_{x\in\\{0,1\\}^{n}}\langle
E_{x}\otimes\tilde{E}_{x}\rvert(P_{k}\otimes\Lambda)\bigg{[}\Big{[}\sum_{k^{\prime}=0}^{2n}\frac{1}{|P_{k^{\prime}}|}\lvert\mathrm{v}(P_{k^{\prime}})\rangle\langle\mathrm{v}(P_{k^{\prime}})\rvert\Big{]}(\mathbb{1}\otimes\Lambda)\Big{[}\sum_{k^{\prime}=0}^{2n}\frac{1}{|P_{k^{\prime}}|}\lvert\mathrm{v}(P_{k^{\prime}})\rangle\langle\mathrm{v}(P_{k^{\prime}})\rvert\Big{]}\bigg{]}^{m}\lvert\rho_{0}\otimes\tilde{\rho}_{0}\rangle.$

Now using the property
$A\otimes\mathbb{1}\lvert\mathrm{v}(B)\rangle=\lvert\mathrm{v}(BA)\rangle$ of
the vectorization function and the orthogonality of the projectors $P_{k}$ we
can rewrite this as

$\displaystyle f_{k}(m)$
$\displaystyle=N_{k}^{-1}\sum_{x\in\\{0,1\\}^{n}}\langle
E_{x}\otimes\Lambda^{\dagger}(\tilde{E}_{x})\rvert\frac{1}{|P_{k^{\prime}}|}\lvert\mathrm{v}(P_{k})\rangle\bigg{[}\frac{1}{|P_{k^{\prime}}|}\langle\mathrm{v}(P_{k})\rvert(\mathbb{1}\otimes\Lambda)\lvert\mathrm{v}(P_{k})\rangle\bigg{]}^{m}\langle\mathrm{v}(P_{k})\rvert\lvert\rho_{0}\otimes\tilde{\rho}_{0}\rangle.$
(75)

Defining

$\displaystyle\lambda_{k}$
$\displaystyle=\frac{1}{|P_{k}|}\operatorname{Tr}(P_{k}\Lambda),$ (76)
$\displaystyle A_{k}$ $\displaystyle=\frac{1}{|P_{k}|}\langle
E_{x}\\!\otimes\\!\Lambda^{\dagger}(\tilde{E}_{x})\lvert\mathrm{v}(P_{k})\rangle\langle\mathrm{v}(P_{k})\rvert\rho_{0}\otimes\tilde{\rho}_{0}\rangle,$
(77)

the expression for $f_{k}(m)$ becomes

$f_{k}(m)=N_{k}^{-1}\sum_{x\in\\{0,1\\}^{n}}A_{k}\lambda^{m}_{k},$ (78)

with $A_{k}=N^{-1}_{k}\sum_{x}A^{x}_{k}$ which is the expression we want.

Let’s now consider the relation between the parameters $\lambda_{k}$ and the
average fidelity. For this is it is useful to remember that the average
fidelity of a trace preserving quantum channel $\Lambda$ can be written as

$F(\Lambda)=\frac{2^{-n}\operatorname{Tr}(\Lambda)+1}{2^{n}+1}.$ (79)

We can explicitly compute this trace from knowledge of the parameters
$\lambda_{k}$, by noting that the projectors $P_{k}$ form a resolution of the
identity (for the space of superoperators), and hence by construction

$\sum_{k=0}^{2n}\frac{|P_{k}|}{2^{2n}}\lambda_{k}=\frac{1}{2^{2n}}\sum_{k=0}^{2n}\operatorname{Tr}(P_{k}\Lambda)=\frac{1}{2^{2n}}\operatorname{Tr}(\Lambda).$
(80)

Inserting this into the above and working out we obtain Eq. 4 in the main
text. Finally we want to consider the parameter $A_{k}$. In order to fit
quantities of the form $A_{k}\lambda_{k}$ it is important that $A_{k}$ is non-
zero. To ensure this we evaluate $A_{k}$ in the noise free limit, setting
$\\{E_{x}\\}$ the be the computational basis and $\rho_{0}$ to be the all-zero
state. The parameter $A_{k}$ is composed of quantities of the form
$\operatorname{Tr}(\lvert x\rangle\\!\langle x\rvert P_{k}(\lvert
x\rangle\\!\langle x\rvert))$. We begin by noting that

$\operatorname{Tr}(\lvert x\rangle\\!\langle x\rvert P_{k}(\lvert
x\rangle\\!\langle x\rvert))=2^{-2n}\sum_{\begin{subarray}{c}S\subset[2n]\\\
|S|=k\end{subarray}}|\operatorname{Tr}(\lvert x\rangle\\!\langle
x\rvert\gamma[S])|^{2}=2^{-2n}\sum_{\begin{subarray}{c}S\subset[2n]\\\
|S|=k\end{subarray}}|\operatorname{Tr}(\lvert 0\rangle\\!\langle
0\rvert\gamma[S])|^{2},$ (81)

since $\langle x\rvert=\langle 0\rvert X_{x}$ for a bit-flip Pauli operator
$X_{x}$ and $X_{x}\gamma[S]X_{x}^{\dagger}=\pm\gamma[S]$. Note that
$|\operatorname{Tr}(\lvert 0\rangle\\!\langle 0\rvert\gamma[S])|$ is one when
$\gamma[S]$ is an all-$Z$ Pauli operator and zero otherwise. From the
definition of the Majorana operators one can see that
$|\operatorname{Tr}(\lvert 0\rangle\\!\langle 0\rvert\gamma[S])|$ is always
zero if $k$ is odd and that for even $k$ there are $\binom{n}{k/2}$ choices
for $S$ such that $\gamma[S]$ is an all-$Z$ Pauli. With all the above we can
evaluate the $A_{k}$ parameters as

$\displaystyle A_{k}$
$\displaystyle=N_{k}^{-1}\sum_{x\in\\{0,1\\}^{n}}\operatorname{Tr}(\lvert
x\rangle\\!\langle x\rvert P_{k}(\lvert x\rangle\\!\langle
x\rvert))\operatorname{Tr}(\lvert 0\rangle\\!\langle 0\rvert P_{k}(\lvert
0\rangle\\!\langle 0\rvert))$ (82)
$\displaystyle=N_{k}^{-1}2^{-n}\frac{\binom{n}{k/2}^{2}}{\binom{2n}{k}},$ (83)

for $k$ even. Hence if we want $A_{k}=1$ in the noise free limit the
normalization in the even $k$, $Z$-basis SPAM case must be set to

$N_{k}^{(Z)}=2^{-n}\frac{\binom{n}{k/2}^{2}}{\binom{2n}{k}},$ (84)

where we also used that $|P_{k}|=\binom{2n}{k}$. We can repeat this exercise
for odd $k$ and $X$-basis SPAM. Here the key expression is

$\displaystyle\operatorname{Tr}\big{(}H^{\otimes n}\lvert x\rangle\\!\langle
x\rvert H^{\otimes n}P_{k}(H^{\otimes n}\lvert x\rangle\\!\langle x\rvert
H^{\otimes n})\big{)}$
$\displaystyle=2^{-n}\sum_{\begin{subarray}{c}S\subset[2n]\\\
|S|=k\end{subarray}}|\operatorname{Tr}\big{(}H^{\otimes n}\lvert
x\rangle\\!\langle x\rvert H^{\otimes n}\gamma[S]\big{)}|^{2}$ (85)
$\displaystyle=2^{-n}\frac{\binom{2n}{(k-1)/2}^{2}}{\binom{2n}{k}},$ (86)

for all $x\in\\{0,1\\}^{n}$. The factor $\binom{2n}{(k-1)/2}$ is obtained by
considering the number of sets $S$ for which
$|\operatorname{Tr}\big{(}H^{\otimes n}\lvert x\rangle\\!\langle x\rvert
H^{\otimes n}\gamma[S]\big{)}|^{2}=1$. This can only happen if $\gamma[S]$
contains only $X$ and $I$ tensor factors. We know that $\gamma_{1}=X_{1}$ and
also that $\gamma_{2i}\gamma_{2i+1}=X_{i}X_{i+1}$ (up to a phase). Since we
require $k$ to be odd the set $S$ must consist of $1$, and $(k-1/2)$ tuples
$2i,2i+1$ with $i\in[1:n]$. Hence there are $\binom{n}{(k-1)/2}$ possible
choices. Using this we can compute

$\displaystyle A_{k}$
$\displaystyle=N_{k}^{-1}\sum_{x\in\\{0,1\\}^{n}}\operatorname{Tr}(H^{\otimes
n}\lvert x\rangle\\!\langle x\rvert H^{\otimes n}P_{k}(H^{\otimes n}\lvert
x\rangle\\!\langle x\rvert H^{\otimes
n}))\operatorname{Tr}(\lvert+\rangle\\!\langle+\rvert
P_{k}(\lvert+\rangle\\!\langle+\rvert))$ (87)
$\displaystyle=N_{k}^{-1}2^{-n}\frac{\binom{n}{\lfloor
k/2\rfloor}^{2}}{\binom{2n}{k}}.$ (88)

Again imposing the unit condition gives the required normalization.

##  Appendix D Gate-dependent noise

In this section we discuss what happens when the gate-independent noise
assumption made above breaks down. While a detailed calculation is out of the
scope of this paper, we aim to show here that the conclusions reached for
randomized benchmarking with arbitrary finite groups [10] essentially carry
over to the case of the matchgate group (which is not finite, but is compact).
Consider a general map $\phi:\mathcal{M}_{n}\to\mathcal{S}_{d}$ from the
orthogonal matchgates to the space of superoperators which assigns to each
orthogonal matchgate a quantum channel. This is a general model of gate-
dependent noise. Note that this is not the most general possible model as we
are ignoring non-Markovian and time-dependent effects. Within this model we
can write the output of the matchgate benchmarking protocol as

$\displaystyle f_{k}(m)$
$\displaystyle=\int_{\operatorname{O}(2n)}dQ_{1}\cdots
Q_{m}\sum_{x\in\\{0,1\\}^{n}}\alpha_{k}(x,Q_{1}\cdots
Q_{m})\langle\tilde{E}_{x}\rvert\phi(Q_{m})\cdots\phi(Q_{1})\lvert\tilde{\rho}_{0}\rangle$
(89) $\displaystyle=N_{k}^{-1}\int_{\operatorname{O}(2n)}dQ_{1}\cdots
Q_{m}\sum_{x\in\\{0,1\\}^{n}}\langle
E_{x}\otimes\tilde{E}_{x}\rvert(P_{k}\otimes\mathbb{1})\omega(Q_{m})\otimes\phi(Q_{m}))\cdots\omega(Q_{1})\otimes\phi(Q_{1})\lvert\rho_{0}\otimes\tilde{\rho}_{0}\rangle,$
(90)

with $\tilde{E}_{x},\tilde{\rho}_{0}$ noisy versions of their ideal
counterparts. Noting that $P_{k}$ commutes with $\omega(Q)$ and that
$P_{k}\omega(Q)=\omega_{k}(Q)$ where $\omega_{k}$ is the relevant irreducible
subrepresentation, we have

$f_{k}(m)=\sum_{x\in\\{0,1\\}^{n}}\langle
E_{x}\otimes\tilde{E}_{x}\rvert\left(\int_{\operatorname{O}(2n)}dQ\,\omega_{k}(Q)\otimes\phi(Q)\right)^{m}\lvert\rho_{0}\otimes\tilde{\rho}_{0}\rangle.$
(91)

We can now consider the operator
$\int_{\operatorname{O}(2n)}dQ\,\omega_{k}(Q)\otimes\phi(Q)$ as a perturbation
of a rank one projector. In that case $f_{k}(m)$ will be well-described by the
$m$-fold power of the largest eigenvalue of
$\int_{\operatorname{O}(2n)}dQ\,\omega_{k}(Q)\otimes\phi(Q)$. A candidate
rank-one projector is given by the ideal implementation
$\int_{\operatorname{O}(2n)}dQ\omega_{k}(Q)\otimes\phi(Q)$. By an application
of Schur’s lemma it can be seen that this operator is equal to $\lvert
v(P_{k})\rangle\langle v(P_{k})\rvert$, a fact we have used already in the
gate-independent derivation. More concretely, assume for a matrix norm
$\|\cdot\|$ on the space $L(\Gamma_{k})\otimes\mathcal{S}_{d}$ that

$\left\|\int_{\operatorname{O}(2n)}dQ\,\omega_{k}(Q)\otimes\phi(Q)-\int_{\operatorname{O}(2n)}dQ\,\omega_{k}(Q)\otimes\omega(Q)\right\|\leq\delta,$
(92)

for some $\delta>0$. From the perturbation theory of $m$-fold matrix powers it
can then be concluded that (provided $\delta$ is small enough)

$f_{k}(m)=A_{k}\lambda_{k}^{m}+O(\delta^{m}),$ (93)

where $\lambda_{k}$ is the largest eigenvalue of the operator
$\int_{\operatorname{O}(2n)}dQ\,\omega_{k}(Q)\otimes\phi(Q)$. In practice this
means that even moderate deviations from the gate-independence assumption get
suppressed exponentially quickly in the sequence length $m$. Hence even when
the gate-independent noise assumption is relaxed the data obtained from a
matchgate benchmarking experiment will be well described by a single
exponential decay.

From the argument above it is not clear how small $\delta$ must be chosen, and
what a physically reasonable choice of submultiplicative norm is. In [39]
(following [10]) it was shown (as a straightforward consequence of their
theorem 1) that Eq. 93 holds for standard randomized benchmarking with a
compact group $\mathbb{G}$ provided

$\int dg\;\|{\phi(g)-\omega(g)}\|_{\diamond}\leq\delta\leq 1/9,$ (94)

holds (where $\|\cdot\|_{\diamond}$ is the diamond norm) and the integral is
taken over the Haar measure. It is important to note here that the factor
$1/9$ is likely suboptimal.

We can get a crude, rule-of-thumb indication for the size of $\delta$ by
considering the decomposition of general matchgates into two-qubit matchgates
given in Eq. 6. Assuming that the single qubit $Z$ rotations are noiseless and
that the two qubit rotations have an average diamond error of $\Delta$ we see
by the triangle inequality and the sub-multiplicativity of the diamond norm
that

$\int dQ\;\|{\phi(Q)-\omega(Q)}\|_{\diamond}\leq n(n-1)\Delta.$ (95)

We can further estimate this by assuming that $\Delta\approx 1-F_{\rm avg}$
where $F_{\rm avg}$ is the average fidelity of the average two qubit gate.
This is of course not generally true, and more or less corresponds to a
"decoherent noise" assumption. In [23] a median two qubit gate fidelity of
$\approx 97\%$ was reported, which we can slot in to give $\delta\approx
n(n-1)0.03$. Hence $\delta\leq 1/9$ for two qubits, and $\delta\approx 0.6$
for five qubits. We emphasize that this is a very crude order of magnitude
estimation (on top of a suboptimal perturbation bound) meant to justify that
$\delta$ can be small compared to $\lambda_{k}$ in reasonable circumstances,
and should not be seen as an upper bound on the tolerance to gate-dependent
noise of our protocol (especially in the context of larger $n$, where the
triangle inequality used above becomes quite wasteful).

##  Appendix E Computation of correlation functions $\alpha_{k}$

We show explicitly how to compute relevant quantities efficiently (in $n$). In
particular the correlation functions given in Eq. 3. For this we will use some
computational techniques from free (or Gaussian) fermionic states and
operations (with which matchgates coincide). For an overview of these
techniques see [20, 40]. We begin by reviewing some identities. For a Majorana
operator $\gamma[S]$ with $S\in[2n],|S|=k$ and a generalized matchgate $U(Q)$,
with $Q\in\operatorname{O}(2n)$ we have

$U(Q)\gamma[S]U(Q)^{\dagger}=\sum_{S^{\prime}\subset[2n],|S^{\prime}|=k}\det\big{(}Q[S,S^{\prime}]\big{)}\gamma[S^{\prime}],$
(96)

where $Q[S,S^{\prime}]$ denotes the matrix $Q$ with only the row indices in
$S$ and column indices in $S^{\prime}$ retained. Moreover we have for a
computational basis state $\lvert x\rangle$ and Majorana $\gamma[S]$ that

$\langle x\rvert\gamma[S]\lvert x\rangle=\operatorname{Pf}(iM_{x}[S]),$ (97)

with $M[S]:=M[S,S]$, and where $\operatorname{Pf}$ denotes the Pfaffian. The
matrix $M_{x}$ is defined by

$M_{x}=\bigoplus_{i=1}^{n}\begin{pmatrix}0&(-1)^{x}_{i}\\\
(-1)^{x_{i}+1}&0\end{pmatrix}.$ (98)

Equation 97 is essentially Wick’s theorem, and it extends to more general
states. The Pfaffian has the following three useful properties:

$\operatorname{Pf}(A)^{2}=\sqrt{\det(A)},$ (99)

for any even dimensional anti-symmetric matrix $A$,

$\operatorname{Pf}\begin{pmatrix}M&C\\\
-C^{T}&N\end{pmatrix}=\operatorname{Pf}(M)\operatorname{Pf}(N+CM^{-1}C^{T}),$
(100)

for invertible $M$, and

$\operatorname{Pf}(QMQ^{T})=\det(Q)\operatorname{Pf}(M).$ (101)

Moreover we will use a more advanced summation identity for Pfaffians, proven
in [24]. Given anti-symmetric $2n\times 2n$ matrices $A,B$ and a $2n\times 2n$
matrix $C$ we have the polynomial identity

$\operatorname{Pf}(A)\operatorname{Pf}\big{(}\frac{A^{\mathrm{co}}}{\operatorname{Pf}(A)}+z^{2}(CBC^{T})\big{)}=\sum_{s=0,s\text{
even}}^{n}z^{s}\sum_{\begin{subarray}{c}S,S^{\prime}\subset[2n]\\\
|S|=|S^{\prime}|=s\end{subarray}}\operatorname{Pf}(A[S])\operatorname{Pf}(B[S^{\prime}])\det(C[S,S^{\prime}]),$
(102)

where ${A}^{\mathrm{co}}$ denotes the co-Pfaffian matrix of $A$, which for
invertible $A$ is given as ${A}^{\mathrm{co}}=\operatorname{Pf}(A)A^{-T}$.
Note that this implies that ${R}^{\mathrm{co}}=\operatorname{Pf}(R)R$ if $R$
is orthogonal. With these identities in hand we move on to compute
$\alpha_{k}(x,Q)$ for even $k$ and $Z$-basis SPAM. From the definition we have

$\displaystyle\alpha_{k}(x,Q)=N_{k}^{-1}\sum_{\vphantom{{{\sum^{i}}^{I}}^{I}}\mathclap{S\subseteq[2n],\lvert
S\rvert=k}}\beta_{S}(x,I)\beta_{S}(0,Q),$

with
$\beta_{S}(x,Q)=2^{-n/2}\operatorname{Tr}\bigl{(}\gamma[S]U(Q)E_{x}U(Q)^{\dagger}\bigr{)}$.
Using $E_{x}=\lvert x\rangle\\!\langle x\rvert$ we see that

$\displaystyle\alpha_{k}(x,Q)=\frac{2^{-n}}{N_{k}}\;\;\;\sum_{\vphantom{{{\sum^{i}}^{I}}^{I}}\mathclap{S\subseteq[2n],\lvert
S\rvert=k}}\operatorname{Pf}(iM_{x}[S])\sum_{\vphantom{{{\sum^{i}}^{I}}^{I}}\mathclap{S^{\prime}\subseteq[2n],\lvert
S^{\prime}\rvert=k}}\operatorname{Pf}(iM_{0}[S^{\prime}])\det(Q[S^{\prime},S]).$

To this we can apply the summation formula Eq. 102 to conclude that

$\alpha_{k}(x,Q)=-\frac{2^{-n}}{N_{k}}\frac{1}{k!}\,\partial^{k}_{z}\left(\operatorname{Pf}(iM_{x})\operatorname{Pf}(iM_{x}+z^{2}iQM_{0}Q^{T})\right)\bigg{|}_{z=0},$
(103)

which can easily be evaluated numerically. Next we calculate the correlator
for odd $k$ and $X$-basis SPAM, which is more complicated. We begin by
calculating the quantity

$\beta_{x}^{Q}[S]=\operatorname{Tr}\big{(}H^{\otimes n}\lvert
x\rangle\\!\langle x\rvert H^{\otimes n}U(Q)\gamma[S]U(Q)^{\dagger}\big{)},$

for $x\in\\{0,1\\}^{n}$ and arbitrary $U(Q)\in\mathcal{M}_{n}^{+}$. We first
note that for every $X$ basis state $H^{\otimes n}\lvert x\rangle$ there
exists $Q_{x}\in\operatorname{O}(2n)$ s.t. $U(Q_{x})H^{\otimes n}\lvert
x\rangle=\frac{I+i\gamma_{1}}{\sqrt{2}}\lvert 0\rangle$. Using this and the
determinant expression for the action of $U(Q)$ we get

$\beta_{x}^{Q}[S]=2^{-1}\;\;\;\sum_{\mathclap{S^{\prime}\subseteq[2n],\lvert
S^{\prime}\rvert=k}}\det\big{(}Q_{x}Q[S,S^{\prime}]\big{)}\lvert
0\rangle(I+i\gamma_{1})\gamma[S^{\prime}](I+i\gamma_{1})\lvert
0\rangle=\sum_{\vphantom{{{\sum^{i}}^{I}}^{I}}\mathclap{S^{\prime}\subseteq[2n],\lvert
S^{\prime}\rvert=k}}\det\big{(}Q_{x}Q[S,S^{\prime}]\big{)}\mathcal{I}(1\in
S^{\prime})\operatorname{Pf}(iM_{0}[S^{\prime}/\\{1\\}]).$

This gives for the correlation function

$\alpha_{k}(x,Q)=\frac{2^{-n}}{N_{k}}\sum_{\begin{subarray}{c}\vphantom{{{\sum^{i}}^{I}}^{I}}\mathclap{S\subseteq[2n],\lvert
S\rvert=k}\\\ \mathclap{S^{\prime}\subseteq[2n],\lvert S^{\prime}\rvert=k}\\\
\mathclap{S^{\prime\prime}\subseteq[2n],\lvert
S^{\prime\prime}\rvert=k}\end{subarray}}\det\big{(}Q_{x}[S,S^{\prime}]\big{)}\det\big{(}Q_{0}Q[S,S^{\prime}]\big{)}\mathcal{I}(1\in
S^{\prime})\mathcal{I}(1\in
S^{\prime\prime})\operatorname{Pf}(iM_{0}[S^{\prime}/\\{1\\}])\operatorname{Pf}(iM_{0}[S^{\prime\prime}/\\{1\\}]),$
(104)

which we can simplify using the Cauchy-Binet identity to

$\alpha_{k}(x,Q)=\frac{2^{-n}}{N_{k}}\sum_{\begin{subarray}{c}\vphantom{{{\sum^{i}}^{I}}^{I}}\mathclap{S^{\prime}\subseteq[2n],\lvert
S^{\prime}\rvert=k}\\\ \mathclap{S^{\prime\prime}\subseteq[2n],\lvert
S^{\prime\prime}\rvert=k}\end{subarray}}\det\big{(}Q_{x}^{T}Q_{0}Q[S^{\prime\prime},S^{\prime}]\big{)}\mathcal{I}(1\in
S^{\prime})\mathcal{I}(1\in
S^{\prime\prime})\operatorname{Pf}(iM_{0}[S^{\prime}/\\{1\\}])\operatorname{Pf}(iM_{0}[S^{\prime\prime}/\\{1\\}]).$
(105)

Now note that $\mathcal{I}(1\in
S^{\prime})\operatorname{Pf}(iM_{0}[S^{\prime}/\\{1\\}])$ is always zero if
$2\in S^{\prime}$ (this follows directly from the definition of $M_{0}$ and
the fact that the Pfaffian is always zero for non-full rank matrices). Hence
we can rewrite the correlation function as

$\alpha_{k}(x,Q)=\frac{2^{-n}}{N_{k}}\sum_{\begin{subarray}{c}\vphantom{{{\sum^{i}}^{I}}^{I}}\mathclap{S^{\prime}\subseteq[3:2n],\lvert
S^{\prime}\rvert=k-1}\\\ \mathclap{S^{\prime\prime}\subseteq[3:2n],\lvert
S^{\prime\prime}\rvert=k-1}\end{subarray}}\det\big{(}Q_{x}^{T}Q_{0}Q[S^{\prime\prime},S^{\prime}]\big{)}\operatorname{Pf}(iM_{0}[S^{\prime}])\operatorname{Pf}(iM_{0}[S^{\prime\prime}]).$
(106)

Defining the matrices
$\widetilde{Q_{x}^{T}Q_{0}Q}=(Q_{x}^{T}Q_{0}Q)\big{[}[3\\!\\!:\\!\\!2n]\big{]}$
and $\widetilde{M_{0}}=M_{0}\big{[}[3\\!\\!:\\!\\!2n]\big{]}$ we can apply Eq.
102 and obtain the correlation function as a $k-1$’th derivative of a Pfaffian
generating function involving matrices of dimension $2(n-1)$.

$\alpha_{k}(x,Q)=\frac{2^{-n}}{N_{k}}\frac{1}{(k-1)!}\partial_{z}^{k-1}\operatorname{Pf}\big{(}i\widetilde{M_{0}}+z^{2}i\widetilde{Q_{x}^{T}Q_{0}Q}\widetilde{M_{0}}(\widetilde{Q_{x}^{T}Q_{0}Q})^{T}\big{)}\bigg{|}_{z=0},$
(107)

which allows for direct numerical calculation.