$|\langle\Psi(d),\Psi\rangle|=1$.

######  Lemma 3.12.

The algorithm $\mathsf{describe}:\mathsf{Mat}_{2n\times
2n}(\mathbb{R})\rightarrow\mathsf{Desc}_{n}$ runs in time $O(n^{3})$. Its
output is such that for every covariance matrix $\Gamma$, the state
$\Psi(\mathsf{describe}(\Gamma))$ is a Gaussian state with covariance matrix
$\Gamma$. We have

$\displaystyle|r|^{2}=|\langle x,\Psi(d)\rangle|^{2}\geq 2^{-n}$ (173)

for $d=(\Gamma,x,r)=\Psi(\mathsf{describe}(\Gamma))$.

###### Proof.

Let $\Gamma\in\mathsf{Mat}_{2n\times 2n}(\mathbb{R})$ be a covariance matrix
and let $\Psi$ be a Gaussian state with covariance matrix $\Gamma$. By
definition of the algorithm $\mathsf{findsupport}$, the value
$y\in\\{0,1\\}^{n}$ computed in line 4 satisfies

$\displaystyle|\langle y,\Psi\rangle|^{2}$ $\displaystyle\geq 2^{-n}\ .$ (174)

By Eq. (45), the value $r$ computed in Line 6 satisfies

$\displaystyle r$ $\displaystyle=|\langle y,\Psi\rangle|\ .$ (175)

In particular, there is an angle $\vartheta\in[0,2\pi)$ such that

$\displaystyle r$ $\displaystyle=\sqrt{\sigma
2^{-n}\mathsf{Pf}(\Gamma(\ket{y})+\Gamma)}=\langle
y,e^{i\vartheta}\Psi\rangle\ .$ (176)

It follows immediately that $d=(\Gamma,y,r)$ is a valid description of the
Gaussian state $e^{i\vartheta}\Psi=\Psi(d)$ with the required property (173).
∎

##  4 Classical simulation of fermionic Gaussian circuits with non-Gaussian
initial states

In this section, we argue that the techniques developed in Section 3 to
describe fermionic Gaussian states (including relative phases) give rise to
efficient classical simulation algorithms for computations composed of non-
Gaussian initial states, Gaussian unitaries and occupation number
measurements. Specifically, we argue that algorithms developed in the context
of stabilizer circuits can immediately be translated to this fermionic setup.
Furthermore, this translation maintains runtime bounds when the stabilizer
extent is replaced by the fermionic Gaussian extent. Because of the generality
of this adaptation procedure – it being applicable to a variety of simulation
algorithms both for strong and weak simulation – we restrict our attention to
the key substitutions.

Our algorithms apply to the efficient classical simulation of fermionic
circuits of the following form, involving $n$ fermions.

  1. (i)

The initial state $\Psi^{(0)}=\Psi$ is a possibly non-Gaussian state $\Psi$.
We assume that its fermionic Gaussian extent $\xi(\Psi)$ and a corresponding
optimal decomposition into a superposition of Gaussian states is known. This
is the case for example for any four-fermion state, or a tensor product of two
four-fermion states, see Section 5. Alternatively, we may assume that an upper
bound $\overline{\xi}(\Psi)\geq\xi(\Psi)$ and a corresponding decomposition of
$\Psi$ achieving this value is known: In this case, runtime upper bounds will
depend on $\overline{\xi}(\Psi)$ instead of $\xi(\Psi)$.

  2. (ii)

The computation proceeds in a sequence of timesteps. At each step $t\in[T]$,
one of the following is performed:

    1. (a)

A Gaussian unitary $U_{R}$, $R\in\mathsf{Gen}(O(2n))$ is applied to the state.
Here the choice of $R$ may depend (in an efficiently computable manner) on
measurement results obtained previously. We will leave this dependence
implicit and do not take it into account in our runtime estimates, as it will
typically depend heavily on the circuit considered.

    2. (b)

An occupation number measurement, i.e., measurement of the operator
$a_{j}^{\dagger}a_{j}$ for some $j\in[n]$ is performed, yielding a measurement
outcome $s\in\\{0,1\\}$ and a corresponding post-measurement state. The choice
of the mode $j\in[n]$ to be measured may again depend (in an efficient manner)
on the measurement outcomes already obtained.

We note that the restriction to the set of Gaussian unitaries associated with
generators of $O(2n)$ in (iia) incurs no loss of generality at the cost of
possibly increasing $T$ by a factor of order $O(n^{2})$ and an additive term
in the runtime of order $O(n^{3})$ since a decomposition of an arbitrary
element $R\in O(2n)$ of the form (36) as a product of $L\leq O(n^{2})$
generators can be found in time $O(n^{3})$, see the discussion below Theorem
3.2.

The use of arbitrary initial states $\Psi$ in (i) allows us to model, in
particular, the application of certain “magic gates” using so-called gadgets.
These can be realized by using non-Gaussian auxiliary states combined with
Gaussian operations, see e.g., [12, 13]. Since all $1$-, $2$\- and $3$-fermion
states are Gaussian [32], $4$-fermion states provide the smallest non-trivial
examples; these will also be our main focus in Section 5. We refer to e.g.,
[12, 13] for a discussion of these constructions.

We proceed as follows: In Section 4, we formulate in general terms how
simulation algorithms for a model can be generalized to initial states that
are superpositions: This follows known approaches for stabilizer circuits
augmented by magic states. In Section 4.2 we review the relationship between
the $\mathcal{D}$-extent and the $\mathcal{D}$-rank defined by a dictionary
$\mathcal{D}$. In Section 4.3 we discuss fast algorithms for estimating norms
of superpositions of dictionary states. In Section 4.4 we apply these
constructions to our setup.

###  4.1 Extending simulation algorithms to superpositions

Here we discuss how to extend simulation algorithms for an efficiently
simulable model $(\mathcal{D},\mathcal{E},\mathcal{M})$ in such a way that the
resulting extended algorithms
$(\chi\mathsf{evolve},\chi\mathsf{measureprob},\allowbreak\chi\mathsf{postmeasure})$
work with any initial state $\Psi$ which is a superposition of $\chi$ elements
of $\mathcal{D}$ (i.e., has $\mathcal{D}$-rank bounded by
$\chi_{\mathcal{D}}(\Psi)\leq\chi$). Our discussion is standard and is
included only for the reader’s convenience: It follows that for stabilizer
states as discussed in [18].

Recall that the dictionary $\mathcal{D}$ is a set of states, $\mathcal{E}$ a
set of operations and $\mathcal{M}$ a set of measurements. In addition to the
subroutines $\mathsf{evolve},\mathsf{measureprob}$ and $\mathsf{postmeasure}$
for evolution and measurement associated with
$(\mathcal{D},\mathcal{E},\mathcal{M})$, the construction discussed here
requires an efficient algorithm $\mathsf{overlap}$ which computes inner
products $\langle\Psi(d_{1}),\Psi(d_{2})\rangle$ from descriptions
$(d_{1},d_{2})\in\mathsf{Desc}^{2}_{n}$. This means that the description
$d\in\mathsf{Desc}_{n}$ of a state $\Psi(d)$ must include phase information.
For Gaussian states, the covariance matrix formalism has to be extended as
discussed in Section 3.

Our goal is to find classical simulation algorithms for circuits of the
following form:

  1. (i)

The initial state $\Psi^{(0)}=\Psi$ is a superposition of the form

$\displaystyle\Psi=\sum_{j=1}^{\chi}\gamma_{j}\varphi_{j}\ .$ (177)

of $\chi$ states $\\{\varphi_{j}\\}_{j=1}^{\chi}\subset\mathcal{D}$ with
complex coefficients $\\{\gamma_{j}\\}_{j=1}^{\chi}$. We assume that this
decomposition is explicitly provided as an input to the classical algorithm in
the form of a $\chi$-tuple $\\{(\gamma_{j},d_{j})\\}_{j=1}^{\mathcal{D}}$,
where $d_{j}$ is an efficient classical descriptions of the state
$\varphi_{j}$.

  2. (ii)

In each timestep $t\in[T]$,

    1. (a)

either an evolution operation $E\in\mathcal{E}$, or

    2. (b)

a measurement $M\in\mathcal{M}$

is applied to the state. We assume that corresponding efficient descriptions
of $E$ respectively $M$ are given to the classical simulation algorithm.

The algorithms
$(\mathsf{evolve},\mathsf{measureprob},\mathsf{postmeasure},\mathsf{overlap})$
associated with the model $(\mathcal{D},\mathcal{E},\mathcal{M})$ then
immediately give rise to algorithms
$(\chi\mathsf{evolve},\chi\mathsf{measureprob},\chi\mathsf{postmeasure})$ for
simulating a more general circuit: At each time step $t\in[T]$, the resulting
algorithm maintains the data
$\\{\gamma_{j}^{(t)},d_{j}^{(t))}\\}_{j=1}^{\mathcal{D}}$ describing the
instantaneous state $\Psi^{(t)}$ after step $t$ as a linear combination

$\displaystyle\Psi^{(t)}$
$\displaystyle=\sum_{j=1}^{\chi}\gamma_{j}^{(t)}\Psi(d_{j}^{(t)})\ $ (178)

of vectors belonging to the dictionary $\mathcal{D}$, and the subroutines
$(\chi\mathsf{evolve},\chi\mathsf{measureprob},\allowbreak\chi\mathsf{postmeasure})$
are used to successively update this description (respectively sample from
corresponding measurement outcomes).

Before describing the extended routines
$\chi\mathsf{evolve},\chi\mathsf{measureprob},\chi\mathsf{postmeasure}$ in
more detail, it is convenient to introduce a subroutine $\chi\mathsf{norm}$
which takes as input a tuple
$\\{(\gamma_{j},d_{j})\\}_{j=1}^{\chi}\in(\mathbb{C}\times\mathsf{Desc}_{n})^{\chi}$
and outputs the squared norm $\|\sum_{j=1}^{\chi}\gamma_{j}\Psi(d_{j})\|^{2}$.
It is clear that such a primitive can be realized naively by using the
algorithm $\mathsf{overlap}$ for computing inner products. This naive
implementation, which we refer to as $\chi\mathsf{naivenorm}$, requires time

$\displaystyle\mathsf{time}(\chi\mathsf{naivenorm})$
$\displaystyle=\chi^{2}\mathsf{time}(\mathsf{overlap})\ .$ (179)

Let us now describe the procedures
$\chi\mathsf{evolve},\chi\mathsf{measureprob}$ and $\chi\mathsf{postmeasure}$,
building on a (general) norm computation subroutine $\chi\mathsf{norm}$.

  1. (a)

if an evolution operation $E\in E$ with description $d_{E}$ is applied at time
$t$, then the description is updated by setting

$\displaystyle\gamma_{j}^{(t)}$
$\displaystyle=\gamma_{j}^{(t-1)}\qquad\textrm{ and }\qquad
d_{j}^{(t)}=\mathsf{evolve}(d_{E},d_{j}^{(t-1)})\qquad\textrm{ for }\qquad
j\in[\chi]\ .$ (180)

This defines the algorithm $\chi\mathsf{evolve}$. The runtime of this update
is

$\displaystyle\mathsf{time}(\chi\mathsf{evolve})=\chi\cdot\mathsf{time}(\mathsf{evolve})\
.$ (181)

  2. (b)

if a (projective) measurement $M=\\{M_{s}\\}_{s\in\mathcal{M}}\in\mathcal{M}$
with description $d_{M}$ is applied to the state at time $t$, then the
probability of obtaining $s\in\mathcal{M}$ is given by

$\displaystyle p(s|\Psi^{(t-1)})$
$\displaystyle=\|M_{s}\Psi^{(t-1)}\|^{2}=\left\|\sum_{j=1}^{\chi}\gamma_{j}^{(t-1)}\sqrt{p(s|\Psi^{(t-1)}_{j})}\Psi^{(t-1)}_{j}(M,s)\right\|^{2}\
.$ (182)

Here the probability
$p(s|\Psi^{(t-1)}_{j})=\|M_{s}\Psi_{j}^{(t-1)}\|^{2}=\mathsf{measureprob}(d^{(t-1)}_{j},d_{M},s)$
of obtaining outcome $s$ when measuring $\Psi^{(t-1)}_{j}$ can be efficiently
obtained from the description $d^{(t-1)}_{j}$ of $\Psi_{j}^{(t-1)}$ and the
description $d_{M}$ of $M$. (Summands $j$ where the probability
$p(s|\Psi_{j}^{(t-1)})$ vanishes can be omitted from this sum.) Similarly, a
description $d_{j}(s)=\mathsf{postmeasure}(d^{(t-1)}_{j},d_{M},s)$ of the
(normalized) post-measurement state
$\Psi^{(t-1)}_{j}(M,s)=\frac{1}{\sqrt{p(s|\Psi^{(t-1)}_{j})}}M_{s}\Psi^{(t-1)}_{j}$
(when measuring $\Psi^{(t-1)}_{j}$) can be obtained efficiently. In
particular, setting
$\tilde{\gamma_{j}}=\gamma^{(t-1)}_{j}\sqrt{p(s|\Psi^{(t-1)}_{j})}$, we
conclude that the outcome probability

$\displaystyle p(s|\Psi^{(t-1)})$
$\displaystyle=\left\|\sum_{j=1}^{\chi}\tilde{\gamma}_{j}\Psi(d_{j}(s))\right\|^{2}\
$ (183)

is the squared norm of a superposition of elements from $\mathcal{D}$. This
expression (together with the norm computation routine $\chi\mathsf{norm}$)
defines the algorithm $\chi\mathsf{measureprob}$. In particular, given
$\\{\tilde{\gamma}_{j},d_{j}(s),p(s|\Psi_{j}^{(t-1)})\\}_{j=1}^{\chi}$, the
probability $p(s|\Psi)$ can be evaluated (exactly) in runtime
$\mathsf{time}(\chi\mathsf{norm})$. Since $\chi\mathsf{measureprob}$ first has
to compute the descriptions $d_{j}(s)$ of the post-measurement states
$\Psi_{j}^{(t-1)}(M,s)$ and the probabilities
$\\{p(s|\Psi_{j}^{(t-1)})\\}_{j=1}^{\chi}$, its runtime is

$\displaystyle\mathsf{time}(\chi\mathsf{measureprob})=\mathsf{time}(\chi\mathsf{norm})+\chi\cdot(\mathsf{time}(\mathsf{measureprob})+\mathsf{time}(\mathsf{postmeasure}))\
.$ (184)

One can easily verify that the post-measurement state after time step $t$ is
given by

$\displaystyle\Psi^{(t)}$
$\displaystyle=\sum_{j=1}^{\chi}\gamma_{j}^{(t)}\Psi(d_{j}^{(t)})\ ,$ (185)

where

$\displaystyle\gamma_{j}^{(t)}=\frac{\tilde{\gamma}_{j}}{\sqrt{p(s|\Psi^{(t-1)})}}\qquad\textrm{
and }\qquad d_{j}^{(t)}=d_{j}(s)\ .$ (186)

In particular, this means that (similarly as for $\chi\mathsf{measureprob}$)
we have an algorithm $\chi\mathsf{postmeasure}$ which given
$\\{(\gamma_{j},d_{j}^{(t-1)})\\}_{j=1}^{\chi}$ and $p(s|\Psi^{(t-1)})$,
computes a description of the post-measurement state in time

$\displaystyle\mathsf{time}(\chi\mathsf{postmeasure})=\chi\cdot(\mathsf{time}(\mathsf{postmeasure})+\mathsf{time}(\mathsf{measureprob}))\
.$ (187)

Given the ability to compute $p(s|\Psi)$ and assuming, e.g., that the number
$|\mathcal{M}|$ of measurement outcomes is constant, one can then sample from
this distribution (when the goal is to perform weak simulation) to get an
outcome $s\in\mathcal{M}$.

Using the naive algorithm $\chi\mathsf{naivenorm}$ for $\chi\mathsf{norm}$
gives runtimes

$\begin{aligned}
\mathsf{time}(\chi\mathsf{evolve})&=\chi\cdot\mathsf{time}(\mathsf{evolve})\\\
\mathsf{time}(\chi\mathsf{measureprob})&=\chi^{2}\cdot\mathsf{time}(\mathsf{overlap})+\chi\cdot(\mathsf{time}(\mathsf{postmeasure})+\mathsf{time}(\mathsf{measureprob}))\\\
\mathsf{time}(\chi\mathsf{postmeasure})&=\chi\cdot(\mathsf{time}(\mathsf{postmeasure})+\mathsf{time}(\mathsf{measureprob}))\end{aligned}\
.$ (188)

As a function of $\chi$, this is dominated by the computation of the squared
norm (183) in $\chi\mathsf{measureprob}$ which takes time $O(\chi^{2})$.

###  4.2 Sparsification: Relating $\mathcal{D}$-extent to approximate
$\mathcal{D}$-rank

Algorithms whose complexity depends on the $\mathcal{D}$-extent
$\xi_{\mathcal{D}}(\Psi)$ instead of the (exact) $\mathcal{D}$-rank
$\chi_{\mathcal{D}}(\Psi)$ (see Eq. (5)) of the initial state $\Psi$ can be
obtained as follows. The idea consists in replacing $\Psi$ by a state
$\tilde{\Psi}$ which is $\delta$-close to $\Psi$ and has bounded
$\mathcal{D}$-rank. More precisely, it relies on the following result which
connects the $\mathcal{D}$-extent $\xi_{\mathcal{D}}(\Psi)$ to the approximate
$\mathcal{D}$-rank $\chi_{\mathcal{D}}^{\delta}(\Psi)$ defined in Eq. (6).

######  Theorem 4.1 (Theorem 1 in [18]).

Suppose $\Psi=\sum_{j=1}^{m}\gamma_{j}\varphi_{j}$ is a decomposition of a
normalized vector $\Psi$ into a superposition of elements
$\\{\varphi_{j}\\}_{j=1}^{m}$ belonging to the dictionary $\mathcal{D}$. Then

$\displaystyle\chi_{\mathcal{D}}^{\delta}(\Psi)\leq
1+\|\gamma\|_{1}^{2}/\delta^{2}$ (189)

where $\|\gamma\|_{1}=\sum_{j=1}^{m}|\gamma_{j}|$ is the $1$-norm of $\gamma$.
In particular, we have the relationship

$\displaystyle\chi_{\mathcal{D}}^{\delta}(\Psi)$ $\displaystyle\leq
1+\xi_{\mathcal{D}}(\Psi)/\delta^{2}\ .$ (190)

In [18], this result was established for the dictionary
$\mathcal{D}=\mathsf{STAB}_{n}$ consisting of $n$-qubit stabilizer states.
Inspection of the proof immediately shows that the statement is applicable to
any dictionary $\mathcal{D}$ (independently of, e.g., whether or not it is
finite). In particular, Theorem 4.1 implies that in runtime upper bounds, the
quantity $\chi_{\mathcal{D}}$ can always be replaced by (the potentially much
smaller quantity) $\xi_{\mathcal{D}}(\Psi)/\delta^{2}$, at the cost of
introducing an $O(\delta)$-error in $L^{1}$-distance in the sampled
distribution. For example, using the naive norm estimation algorithm (i.e.,
inserting into (188)), this gives a quadratic scaling (for computing output
probabilities) in $\xi_{\mathcal{D}}(\Psi)$. Note that here we are assuming
that a decomposition of $\Psi$ with squared $L^{1}$-norm $\|\gamma\|_{1}^{2}$
of coefficients achieving $\xi_{\mathcal{D}}(\Psi)$ is given.

###  4.3 Fast norm estimation and approximate simulation

In pioneering work, Bravyi and Gosset [17] improved upon the sketched
algorithm in the case of stabilizer states. This was achieved by replacing the
$O(\chi^{2})$-runtime (naive) estimation algorithm $\chi\mathsf{naivenorm}$
for the norm of a superposition of stabilizer states by a probabilistic
algorithm $\chi\mathsf{fastnorm}$. With success probability at least
$1-p_{f}$, the algorithm $\chi\mathsf{fastnorm}$ provides an estimate
$\hat{N}$ of the squared norm $N=\left\|\Psi\right\|^{2}$ of a superposition
$\Psi=\sum_{j=1}^{\chi}\gamma_{j}\varphi_{j}$ of $n$-qubit stabilizer states
$\\{\varphi_{j}\\}_{j=1}^{\chi}$ with multiplicative error $\epsilon$ (i.e.,
$\hat{N}\in[(1-\epsilon)N,(1+\epsilon)N])$, and has runtime $O(\chi\cdot
n^{3}\epsilon^{-2}p_{f}^{-1})$ The key observation underlying the algorithm is
the fact that the norm of interest can be expressed as

$\displaystyle N$
$\displaystyle=2^{n}\mathbb{E}_{\Theta}\left[|\langle\Theta,\Psi\rangle|^{2}\right]\
,$ (191)

i.e., it is proportional to the expected (squared) overlap of $\Psi$ with a
state $\ket{\Theta}$ drawn uniformly at random from the set of $n$-qubit
stabilizer states. Given $\\{\gamma_{j},\varphi_{j}\\}_{j=1}^{\chi}$, this
algorithm proceeds by taking $R=\lceil p_{f}^{-1}\epsilon^{-2}\rceil$
stabilizer states $\Theta_{1},\ldots,\Theta_{R}$ chosen uniformly from the set
of all stabilizer states, and producing the estimate

$\displaystyle\hat{N}=\frac{2^{n}}{R}\sum_{k=1}^{R}|\langle\Theta_{k},\Psi\rangle|^{2}\
$ (192)

of $N$. Importantly, expression (192) can be computed from (the descriptions
of) $\\{\Theta_{k}\\}_{k=1}^{R}$, $\\{\varphi_{j}\\}_{j=1}^{\chi}$ and the
coefficients $\\{\gamma_{j}\\}_{j=1}^{\chi}$, using $\chi\cdot R$ calls of a
subroutine $\mathsf{overlap}$ which computes the overlap of two stabilizer
states (including phases). This is because each summand in (192) can be
written as a sum

$\displaystyle|\langle\Theta_{k},\Psi\rangle|$
$\displaystyle=\left|\sum_{j=1}^{\chi}\gamma_{j}\langle\Theta_{k},\varphi_{j}\rangle\right|$
(193)

of $\chi$ such products. The resulting runtime of this norm estimation
algorithm is thus $O(\chi\cdot R\cdot\mathsf{time}(\mathsf{overlap}))$ which
amounts to the claimed runtime $O(\chi\cdot n^{3}\epsilon^{-2}p_{f}^{-1})$.

We note that a similar reasoning can be applied to any situation where the
norm of a superposition of dictionary elements of interest can be expressed as
in Eq. (191) as the expected overlap of the inner product of $\Psi$ with a
state $\Theta$ randomly chosen according to a suitable distribution over
dictionary states. Specifically, as derived in Appendix B of [60] and
discussed below (see Section 4.4), this is the case for the set of fermionic
Gaussian states. The corresponding norm estimation algorithm then has a
runtime of the form

$\displaystyle\mathsf{time}(\chi\mathsf{fastnorm})$ $\displaystyle=O(\chi\cdot
R\cdot\mathsf{time}(\mathsf{overlap})+R\cdot\mathsf{time}(\mathsf{samplestate}))$
(194)
$\displaystyle=O(p_{f}^{-1}\epsilon^{-2}(\chi\mathsf{time}(\mathsf{overlap})+\mathsf{time}(\mathsf{samplestate}))$

where $\mathsf{samplestate}$ is an algorithm producing a description of a
state $\Theta$ drawn randomly form the appropriate distribution. Importantly,
the runtime (194) is linear in $\chi$, resulting in a linear dependence when
replacing $\chi$ by the stabilizer extent $\xi_{\mathcal{D}}(\Psi)$ as
discussed in Section 4.2.

Algorithms
$(\mathsf{approx}\mathsf{evolve},\mathsf{approx}\mathsf{measureprob},\mathsf{approx}\mathsf{postmeasure})$
can now be obtained by using $\chi\mathsf{fastnorm}$ in place of
$\chi\mathsf{norm}$. The algorithm $\mathsf{approx}\mathsf{evolve}$ is
identical to $\chi\mathsf{evolve}$ since it does not involve norm
computations. In contrast, $\mathsf{approx}\mathsf{measureprob}$ is a
probabilistic algorithm that can fail with probability $p_{f}$ and both
$\mathsf{approx}\mathsf{measureprob}$ and
$\mathsf{approx}\mathsf{postmeasure}$ introduce an error (in the sampled
distribution and the post-measurement state, respectively). This is because
$\chi\mathsf{fastnorm}$ only estimates the norm of a superposition.

Finally, replacing $\chi$ by the $\mathcal{D}$-extent
$\xi_{\mathcal{D}}(\Psi)$, see Section 4.2 results in a triple of approximate
algorithms
$(\mathsf{approx}\mathsf{evolve},\mathsf{approx}\mathsf{measureprob},\mathsf{approx}\mathsf{postmeasure})$
with parameters $(\epsilon,\delta,p_{f})$ describing the quality of
approximation and failure probability as discussed in Section 1.5. By
construction, the runtimes of these algorithms are

$\begin{aligned}
\mathsf{time}(\mathsf{approx}\mathsf{evolve})&=O\left(\frac{\xi_{\mathcal{D}}(\Psi)}{\delta^{2}}\mathsf{time}(\mathsf{evolve})\right)\\\
\mathsf{time}(\mathsf{approx}\mathsf{measureprob})&=O\left(p_{f}^{-1}\epsilon^{-2}\left(\frac{\xi_{\mathcal{D}}(\Psi)}{\delta^{2}}\mathsf{time}(\mathsf{overlap})+\mathsf{time}(\mathsf{samplestate})\right)\right)\\\
&+O\left(\frac{\xi_{\mathcal{D}}(\Psi)}{\delta^{2}}\left(\mathsf{time}(\mathsf{postmeasure})+\mathsf{time}(\mathsf{measureprob})\right)\right)\\\
\mathsf{time}(\mathsf{approx}\mathsf{postmeasure})&=O\left(\frac{\xi_{\mathcal{D}}(\Psi)}{\delta^{2}}\left(\mathsf{time}(\mathsf{postmeasure})+\mathsf{time}(\mathsf{samplestate})\right)\right)\end{aligned}\
.$ (195)

###  4.4 Fermionic linear optics with non-Gaussian initial states

The algorithms described above can be adapted in a straightforward manner to
the problem of classically simulating fermionic linear optics with non-
Gaussian initial states: We can simply use the efficient description of
Gaussian states introduced in Section 3 and make use of the associated
procedures $\mathsf{overlap}$, $\mathsf{evolve}$ as well as
$\mathsf{measureprob}$ and $\mathsf{postmeasure}$. In particular, observe that
combining Eq. (188) with the runtimes $O(n^{3})$ for the algorithms
$\mathsf{evolve}$, $\mathsf{postmeasure}$ and $\mathsf{overlap}$, and $O(1)$
for $\mathsf{measureprob}$ (see Section 3) results in the runtimes given in
Table 3 for exact simulation.

To get a linear scaling in the Gaussian extent $\xi_{\mathcal{G}_{n}}(\Psi)$
of the initial state (for approximate simulation), the naive norm estimation
needs to be replaced. A fast norm estimation scheme for superpositions of
fermionic Gaussian has been described in Appendix C of Ref. [60]: Consider the
following probabilistic process defined for a superposition
$\Psi=\sum_{j=1}^{\chi}\gamma_{j}\varphi_{j}$,
$\varphi_{j}\in\mathcal{G}_{n}$, $\gamma_{j}\in\mathbb{C}$ of $n$-mode
fermionic Gaussian states:

  1. (i)

Sample $K$ random Gaussian states $\\{\Theta_{k}\\}_{k=1}^{K}$ independently
and identically from the distribution induced by picking a permutation $\pi\in
S_{2n}$ and a string $y\in\\{0,1\\}^{n}$ uniformly at random and outputting

$\displaystyle\ket{\Theta(\pi,y)}=U_{R_{\pi}}\ket{y}.$ (196)

Here $R_{\pi}=O(2n)$ is a permutation matrix specified by an element $\pi\in
S_{2n}$ and $y\in\\{0,1\\}^{n}$.

  2. (ii)

Set

$\displaystyle\hat{N}=\frac{1}{K}\sum_{k=1}^{K}2^{n}|\langle\Theta_{k},\Psi\rangle|^{2}\
.$ (197)

The following was shown in [60].

######  Lemma 4.2 (Lemma 10 in Ref. [60]).

For any $p_{f}\in[0,1]$ and $\epsilon>0$, consider the probabilistic process
described above with the choice $K=\lceil
2\sqrt{n}\epsilon^{-2}p_{f}^{-1}\rceil$. Then the random variable $\hat{N}$
satisfies

$(1-\epsilon)\|\Psi\|^{2}\leq\hat{N}\leq(1+\epsilon)\|\Psi\|^{2}\ .$

with probability at least $1-p_{f}$.

A description of a state proportional to $\Theta_{k}$ can be computed from the
associated pair $(\pi_{k},y_{k})\in S_{2n}\times\\{0,1\\}^{n}$ using the
subroutine $\mathsf{describe}$, for each $k\in[K]$, as follows (see Fig. 15
for pseudocode for the algorithm). The covariance matrix
$\Gamma_{k}=R_{\pi_{k}}\Gamma(\ket{y})R_{\pi_{k}}^{\dagger}$ of such a state
can be computed in time $O(n^{3})$ from $(\pi_{k},y)$, and applying
$\mathsf{describe}$ to $\Gamma_{k}$ gives the desired description. We note
that any such state can be used in place of $\Theta_{k}$ since the expression
(197) (and, in particular, (193)) does not depend on the global phase of
$\Theta_{k}$. With the definition (197), it follows that the probabilistic
process described here can be simulated in time given by Eq. (194) using $K$
calls to the subroutine $\mathsf{samplestate}$ shown in Fig. 15, and
subsequent use of $\mathsf{overlap}$ to compute the empirical average (197).

1:function $\mathsf{samplestate}$ 2: $y\leftarrow$ uniform random string in
$\\{0,1\\}^{n}$ 3: $\pi\leftarrow$ uniform random permutation in $S_{2n}$ 4:
$R_{\pi}\in O(2n)$ 5: for $j\in[2n]$ do $\triangleright$ compute the
permutation matrix $R_{\pi}$ 6: $R_{\pi}[j]\leftarrow e_{\pi(j)}$ 7:
$\Gamma\leftarrow R_{\pi}\Gamma(\ket{y]})R_{\pi}^{\dagger}$ $\triangleright$
compute the covariance matrix of $\ket{\Theta(\pi,y)}=R_{\pi}\ket{y}$ 8:
return $\mathsf{describe}(\Gamma)$ $\triangleright$ output a description of
$\ket{\Theta(\pi,y)}$ Figure 15: The algorithm $\mathsf{samplestate}$ outputs
a classical description of a state $\ket{\Theta(\pi,y)}=R_{\pi}\ket{y}$ where
$\pi\in S_{2n}$ and $y\in\\{0,1\\}^{n}$ are taken uniformly at random.

Because the runtimes of $\mathsf{describe}$ and $\mathsf{overlap}$ are both
upper bounded by $O(n^{3})$, this leads to an overall runtime of
$O\left(n^{7/2}\epsilon^{-2}p_{f}^{-1}\chi\right)$ of this algorithm for
computing the estimate $\hat{N}$ of $\|\Psi\|^{2}$. We note this conclusion
about the runtime was also reached in Ref. [60], although the issue of a
potential lack of a phase reference applicable throughout the computation was
not considered there. This issue is resolved by our description of Gaussian
states, see Section 3.

Combing this algorithm with the runtimes given in Eq. (195) and with the
runtimes $O(n^{3})$ for the algorithms $\mathsf{evolve}$,
$\mathsf{postmeasure}$ and $\mathsf{overlap}$, and $O(1)$ for
$\mathsf{measureprob}$ (see Section 3) gives runtimes claimed in Table 4 for
the algorithms $\mathsf{approx}\mathsf{evolve}$,
$\mathsf{approx}\mathsf{measureprob}$ and
$\mathsf{approx}\mathsf{postmeasure}$.

To give an idea of the feasibility of running these algorithms, let us compare
to [17], where Clifford circuits on $n=40$ qubits with up to $t=48$ $T$-gates
were simulated by means of an algorithm relying on the stabilizer formalism
augmented with stabilizer-decompositions of magic states. These simulations
were performed using MATLAB on a laptop with a 2.6GHz Intel i5 Dual Core CPU.
Runtimes of individual building blocks are given in [17, Table I] and are of
the order of tens of milliseconds even for 100 qubits. The algorithms of [17]
have a scaling which is linear in the approximate (stabilizer) rank. Since
this quantity scales exponentially in the number of copies of a single magic
state, this dominates the runtime, which includes additional factors that are
polynomial in the number of qubits and gates. In more detail, the authors of
[17] consider the non-stabilizer state

$\displaystyle|A\rangle=\frac{1}{\sqrt{2}}(|0\rangle+e^{i\pi/4}|1\rangle)\ .$
(198)

A copy of this (magic) state can be used together with stabilizer operations
to implement the (non-Clifford) $T$-gate [11]. The approximate stabilizer rank
(which is defined in terms of fidelity rather than distance in [17]) of $t$
copies of this state is shown to be bounded by

$\displaystyle\tilde{\chi}^{\delta}(|A\rangle^{\otimes t})=O(2^{\gamma
t}\delta^{-1})\qquad\textrm{ where }\qquad\gamma\approx 0.228\ .$ (199)

Now consider our simulation procedure: Here the dependence of the runtime is
linear in the Gaussian extent $\xi_{\mathcal{G}_{n}}$ and thus again
essentially linear in the $\delta$-approximate Gaussian rank
$\chi_{\mathcal{G}_{n}}^{\delta}$ (see Theorem 4.1), and this quantity
dominates the runtime (which includes polynomial factors in the number of
fermionic modes and gates). Thus our simulation algorithms have a runtime
comparable to that of [17] when substituting the fermionic Gaussian extent of
the initial state with the approximate stabilizer rank (respectively the
stabilizer extent) of the initial state in the Clifford setting.

To make this more concrete, consider the state

$\displaystyle|a_{8}\rangle=\frac{1}{\sqrt{2}}\left(|0000\rangle+|1111\rangle\right)\
,$ (200)

which is a magic state for Gaussian computations: As shown in [12, Lemma 1],
the non-Gaussian unitary (magic gate) $U_{\textrm{magic}}=\exp(i\pi
c_{1}c_{2}c_{3}c_{4})$ can be implemented with Gaussian operations and a
single copy of this state. The state $|a_{8}\rangle$ has Gaussian extent

$\displaystyle\xi_{\mathcal{G}_{n}}(|a_{8}\rangle)=2\ $ (201)

and thus we have

$\displaystyle\xi_{\mathcal{G}_{n}}(|a_{8}\rangle^{\otimes t^{\prime}})\leq
2^{t^{\prime}}$ (202)

for $t^{\prime}$ copies of $|a_{8}\rangle$. (We note that in fact, the results
of [49] imply that equality holds in (202).) Eq. (201) can be shown as
follows: The upper bound $\xi_{\mathcal{G}_{n}}(|a_{8}\rangle)\leq 2$ is a
direct consequence of the decomposition given in (200) into a sum of the
Gaussian states $|0000\rangle$ and $|1111\rangle$. The dual formulation of the
extent implies that

$\displaystyle\xi_{\mathcal{G}_{n}}(|\Psi\rangle)$
$\displaystyle\geq\frac{1}{F_{\mathcal{G}_{n}}(|\Psi\rangle)}\qquad\textrm{
for any state }\Psi\in\mathcal{H}_{n}\ ,$ (203)

see (294) with the choice
$y=|\Psi\rangle/\sqrt{F_{\mathcal{G}_{n}}(|\Psi\rangle)}$. Using the fact that
$F_{\mathcal{G}_{n}}(|a_{8}\rangle)=1/2$ (see [50] – the inequality
$F_{\mathcal{G}_{n}}(|a_{8}\rangle)\leq 1/2$ which is needed here was
previously shown in [32]) and (203) gives the lower bound
$\xi_{\mathcal{G}_{n}}(|a_{8}\rangle)\geq 2$, establishing (201).

Comparing the exponents in (199) and (202), we can give the following rough
estimate: we expect that the cost (runtime) of simulating $t=48$ $T$-gates in
the stabilizer framework is comparable with the cost of simulation of roughly
$t^{\prime}\approx 10$ gates $U_{\textrm{magic}}$ in the fermionic context.
Similar estimates apply to other non-Gaussian operations such that the SWAP
gate, a magic gate in the context of fermionic Gaussian computation. The
latter can be implemented using a certain non-Gaussian state as shown in [13].

###  4.5 Efficient additive-error strong simulation

In a different direction of generalization, building upon the work [17] and
making innovative use of a concentration inequality by Hayes [62] for vector
martingales, Ref. [27] gives a randomized algorithm which, for a state $\Psi$
obtained by applying $n$-qubit Clifford gates and $t$ (non-Clifford) $T$-gates
to $\ket{0}^{\otimes n}$, provides an additive-error estimate $\hat{p}(x)$ for
the probability $p(x)=\|(\langle x|\otimes
I^{\otimes(n-a)})|\Psi\rangle\|^{2}$ of observing $a$ qubits in the state
$\ket{x}$, with $x\in\\{0,1\\}^{a}$. The algorithm is based on a procedure by
which the probability $p(x)$ of interest is expressed in terms of the squared
norm $\left\|\left(\bra{0}^{\otimes t-r}\otimes I^{\otimes
r}\right)W\ket{\Psi}\right\|^{2}$ of a partially projected state, where $\Psi$
is a tensor product of $t$ non-stabilizer single-qubit states (arising from
gadgetization), $W$ a certain Clifford unitary, and $r$ a circuit-dependent
integer. The failure probability of the constructed algorithm is then upper
bounded (see [27, Theorem 3]) by an expression depending on $p(x)$, the error
$\epsilon$, the stabilizer rank $\xi_{\mathsf{STAB}_{n}}(\Psi)$ (taking the
role of $\chi$) of the product state $\Psi$, as well as two additional
parameters than be chosen freely, but enter into the (polynomial) runtime
estimate.

The described method of adapting fast algorithms for simulating Clifford
circuits with non-stabilizer initial states can be applied in a similar manner
to this algorithm, since this also reduces to computing inner products
(including phases) between Gaussian states.

##  5 Multiplicativity of the Gaussian fidelity for $4$ fermions

The main result of this section is a proof that the fermionic Gaussian
fidelity is multiplicative for the tensor product of any two positive-parity
4-fermion states.

We begin in Section 5.1 by laying out some specific properties of $4$-fermion
states. We discuss a Gaussianity condition specific to 4-fermion states [58]
and we write an explicit expression for any 4-fermion state as a superposition
of two orthogonal (Gaussian) states. This was first introduced in Refs. [58,
56]. In Section 5.2 we establish properties of the fermionic Gaussian fidelity
for 4-fermion states which are subsequently used in Section 5.3 to prove that
the fermionic Gaussian fidelity is multiplicative for the tensor product of
any two 4-fermion states.

###  5.1 Four-fermion Gaussian and non-Gaussian states

Key to our considerations is a certain antiunitary map $\theta$ acting on
$\mathcal{H}_{+}^{4}$, the positive-parity subspace of 4 fermions spanned by
$\\{\ket{x}\\}_{x\in\\{0,1\\}^{4}_{+}}$. It is defined by its action

$\displaystyle\theta\ket{x}$
$\displaystyle=(-1)^{\vartheta(x)}\ket{\overline{x}}\ ,$ (204)

for $x\in\\{0,1\\}_{+}^{4}$, on basis states (antilinearly extended to all of
$\mathcal{H}_{+}^{4}$), where $\vartheta(x)=x_{1}+x_{3}\mod 2=x_{2}+x_{4}\mod
2=\vartheta(\bar{x})$. Here
$\overline{x}=(\overline{x}_{1},\ldots,\overline{x}_{n})$ is obtained by
flipping each bit of $x$. The relevant properties of this map are the
following. We note that the following statement has been given in [58, Eq.
(9)], along with a negative-parity version.

######  Lemma 5.1 ([58]).

A state $\Psi\in\mathcal{H}_{+}^{4}$ is Gaussian if and only if

$\displaystyle\langle\Psi,\theta\Psi\rangle$ $\displaystyle=0\ .$ (205)

###### Proof.

This follows from the Gaussianity criterion given in Lemma 2.3. We give the
proof in Appendix A. ∎

######  Lemma 5.2.

We have $\theta c_{j}c_{k}=c_{j}c_{k}\theta~{}$ for all $j,k\in[8]$.

###### Proof.

See Appendix B. ∎

The following result was first shown in Ref. [56].

######  Lemma 5.3 ([56]).

Let $\Psi\in\mathcal{H}_{+}^{4}$ be a unit vector. Then there are two
orthogonal unit vectors $\Psi_{1},\Psi_{2}\in\mathcal{H}_{+}^{4}$,
$\varphi\in[0,2\pi)$ and $a\in[0,1/\sqrt{2}]$ such that

$\displaystyle\theta\Psi_{j}$ $\displaystyle=\Psi_{j}\qquad\textrm{ for
}\qquad j\in[2]$ (206)

and

$\displaystyle\Psi$
$\displaystyle=e^{i\varphi}\left(\sqrt{1-a^{2}}\Psi_{1}+ia\Psi_{2}\right)\ .$
(207)

###### Proof.

We first argue that it suffices to consider the case where $\Psi$ satisfies

$\displaystyle\langle\Psi,\theta\Psi\rangle\in\mathbb{R}\ .$ (208)

This is because

$\displaystyle\langle(e^{i\varphi}\Psi),\theta(e^{i\varphi}\Psi)\rangle$
$\displaystyle=e^{-2i\varphi}\langle\Psi,\theta\Psi\rangle\qquad\textrm{ for
every }\qquad\varphi\in[0,2\pi)\ ,$ (209)

which implies that (208) can be ensured by replacing $\Psi$ with
$e^{i\varphi}\Psi$ for a suitably chosen $\varphi\in[0,2\pi)$.

Let $\Psi$ be such that (208) holds. We define “real” and “imaginary” parts of
$\Psi$ by the expressions

$\displaystyle\Psi_{R}$
$\displaystyle=\frac{1}{2}\left(\Psi+\theta\Psi\right)$ (210)
$\displaystyle\Psi_{I}$
$\displaystyle=\frac{1}{2i}\left(\Psi-\theta\Psi\right)\ .$ (211)

It follows immediately from this definition that

$\displaystyle\Psi$ $\displaystyle=\Psi_{R}+i\Psi_{I}\ $ (212)

and

$\displaystyle\theta\Psi_{R}$ $\displaystyle=\Psi_{R}$ (213)
$\displaystyle\theta\Psi_{I}$ $\displaystyle=\Psi_{I}\ .$ (214)

because $\theta$ is antiunitary and an involution. Furthermore, Eq. (208)
implies that the vectors $\Psi_{R}$,$\Psi_{I}$ are orthogonal: We have

$\displaystyle 4i\langle\Psi_{R},\Psi_{I}\rangle$
$\displaystyle=\langle\Psi+\theta\Psi,\Psi-\theta\Psi\rangle$ (215)
$\displaystyle=\|\Psi\|^{2}-\|\theta\Psi\|^{2}+\langle\theta\Psi,\Psi\rangle-\langle\Psi,\theta\Psi\rangle$
(216) $\displaystyle=-2i\mathsf{Im}\langle\Psi,\theta\Psi\rangle$ (217)
$\displaystyle=0$ (218)

where we used that $\theta$ is an antiunitary in the first step, and
assumption (208) in the last step. The claim now follows by setting

$\displaystyle\left(a,\Psi_{1},\Psi_{2}\right)=\begin{cases}\left(\|\Psi_{I}\|,\frac{\Psi_{R}}{\|\Psi_{R}\|},\frac{\Psi_{I}}{\|\Psi_{I}\|}\right)\qquad\textrm{
if }\qquad\|\Psi_{I}\|\leq 1/\sqrt{2}\\\
\left(\|\Psi_{R}\|,\frac{\Psi_{I}}{\|\Psi_{I}\|},\frac{\Psi_{R}}{\|\Psi_{R}\|}\right)\qquad\textrm{
otherwise }\ .\end{cases}$ (219)

∎

######  Theorem 5.4 ([58, 56]).

Let $\Psi\in\mathcal{H}_{+}^{4}$ be an arbitrary unit vector. Then there are a
Gaussian pure state $\Psi_{g}\in\mathcal{G}^{+}_{4}$, $\varphi\in[0,2\pi)$ and
$f\in[1/2,1]$ such that the state $\theta\Psi_{g}$ is Gaussian and orthogonal
to $\Psi_{g}$ and

$\displaystyle\Psi$
$\displaystyle=e^{i\varphi}\left(\sqrt{f}\Psi_{g}+\sqrt{1-f}\theta\Psi_{g}\right)\
.$ (220)

The triple $\left(\Psi_{g},\varphi,f\right)$ is uniquely defined by $\Psi$,
i.e., a function of $\Psi$. Furthermore, the quantity $f=f(\Psi)$ is invariant
under the action of Gaussian unitaries associated with special orthogonal
rotations: We have

$\displaystyle f(U\Psi)$ $\displaystyle=f(\Psi)\qquad\textrm{ for any Gaussian
unitary }U=U_{R}\text{ with }R\in SO(2n)\ .$ (221)

###### Proof.

Let $\Psi\in\mathcal{H}_{+}^{4}$ be an arbitrary unit vector. Let
$\Psi_{1},\Psi_{2}\in\mathcal{H}_{+}^{4}$, $\varphi\in[0,2\pi)$ and
$a\in[0,1]$ be as in Lemma 5.3. Define

$\displaystyle\Psi_{g}^{\pm}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\Psi_{1}\pm i\Psi_{2}\right)\ .$ (222)

Then

$\displaystyle\theta\Psi_{g}^{-}=\Psi_{g}^{+}$ (223)

because of property (206). It follows that

$\displaystyle\langle\Psi_{g}^{-},\Psi_{g}^{+}\rangle$
$\displaystyle=\langle\Psi_{g}^{-},\theta\Psi_{g}^{-}\rangle$ (224)
$\displaystyle=\frac{1}{2}\left(\|\Psi_{1}\|^{2}+2i\mathsf{Re}\langle\Psi_{1},\Psi_{2}\rangle-\|\Psi_{2}\|^{2}\right)$
(225) $\displaystyle=0$ (226)

since $\Psi_{1}$ and $\Psi_{2}$ are orthogonal unit vectors.

Using (223) and the orthogonality of $\Psi_{g}^{-},\Psi_{g}^{+}$ implies that

$\displaystyle\langle\Psi_{g}^{+},\theta\Psi_{g}^{+}\rangle$
$\displaystyle=\langle\Psi_{g}^{+},\Psi_{g}^{-}\rangle=0$ (227)

and similarly (because $\theta$ is an involution)

$\displaystyle\langle\Psi_{g}^{-},\theta\Psi_{g}^{-}\rangle$
$\displaystyle=\langle\Psi_{g}^{-},\Psi_{g}^{+}\rangle=0\ .$ (228)

According to the Gaussianity criterion in Lemma 5.1, we conclude that both
$\Psi_{g}^{+}$ and $\Psi_{g}^{-}$ are Gaussian.

Rewriting Eq. (207) by expressing $(\Psi_{1},\Psi_{2})$ in terms of
$(\Psi_{g}^{+},\Psi_{g}^{-})$ gives

$\displaystyle\Psi$
$\displaystyle=e^{i\varphi}\left(\sqrt{f}\Psi_{g}^{-}+\sqrt{1-f}\Psi_{g}^{+}\right)\qquad\textrm{
where }\qquad f=\frac{1}{2}+a\sqrt{1-a^{2}}\ .$ (229)

The claim (220) now follows with (223).

It remains to show property (221) of the function $f$. This follows
immediately from the fact that the antiunitary $\theta$ commutes with all
quadratic monomials $c_{j}c_{k}$ of Majorana operators (see Lemma 5.2), and
hence with any Gaussian unitary $U=U_{R}$ with $R\in SO(2n)$, i.e.,
$U\theta=\theta U$. Retracing the steps of the proof, it is easy to check that
if $(\Psi_{1},\Psi_{2})$ are the states of Lemma 5.3, and $\Psi_{g}$ the state
in expression (220) for $\Psi$, then the corresponding states
$(\Psi^{\prime}_{1},\Psi^{\prime}_{2})$ and $\Psi^{\prime}_{g}$ for the state
$\Psi^{\prime}=U\Psi$ are given by $\Psi^{\prime}_{j}=U\Psi_{j}$ for $j\in[2]$
and $\Psi_{g}^{\prime}=U\Psi_{g}$, respectively. This implies the claim. ∎

###  5.2 The Gaussian fidelity for 4-fermion states

For a subset $E\subset\\{0,1\\}^{4}$, we define
$\overline{E}=\left\\{\overline{x}\ |\ x\in E\right\\}$. We also write
$\Pi_{E}=\sum_{x\in E}|x\rangle\langle x|$ for the projection onto the span of
$\\{\ket{x}\\}_{x\in E}$.

######  Lemma 5.5.

Let $E\subset\\{0,1\\}_{+}^{4}$, $|E|=4$ be a subset of even-weight strings
such that $E\cup\overline{E}=\\{0,1\\}^{4}_{+}$. Let $f(\Psi)\in[1/2,1]$ be
defined as in Theorem 5.4. Then

$\displaystyle\left\|\Pi_{E}\Psi\right\|^{2}\leq f(\Psi)\qquad\textrm{ for any
unit vector }\qquad\Psi\in\mathcal{H}_{+}^{4}\ .$ (230)

###### Proof.

Let $f=f(\Psi)\in[1/2,1]$, $\varphi\in[0,2\pi)$ and
$\Psi_{g}\in\mathcal{G}^{+}_{n}$ be as in Theorem 5.4 such that

$\displaystyle\Psi$
$\displaystyle=e^{i\varphi}\left(\sqrt{f}\Psi_{g}+\sqrt{1-f}\theta\Psi_{g}\right)\
.$ (231)

We define

$\displaystyle\begin{matrix}\alpha_{x}&=&\langle x,\Psi_{g}\rangle\\\
\beta_{x}&=&\langle x,\theta\Psi_{g}\rangle\end{matrix}\qquad\textrm{ for
every }\qquad x\in E\ .$ (232)

We claim that we have the identities

$\displaystyle\sum_{x\in E}(|\alpha_{x}|^{2}+|\beta_{x}|^{2})$
$\displaystyle=1$ (233) $\displaystyle\sum_{x\in
E}\overline{\alpha_{x}}\beta_{x}$ $\displaystyle=0\ .$ (234)

Observe that these two identities immediately imply (230): Using expression
(231), we have

$\displaystyle\|\Pi_{E}\Psi\|^{2}$ $\displaystyle=\sum_{x\in E}|\langle
x,\Psi\rangle|^{2}$ (235) $\displaystyle=\sum_{x\in
E}|\sqrt{f}\alpha_{x}+\sqrt{1-f}\beta_{x}|^{2}$ (236)
$\displaystyle=\left\|\sqrt{f}\vec{\alpha}+\sqrt{1-f}\vec{\beta}\right\|^{2}\
$ (237)

where we defined the vectors $\vec{\alpha}=(\alpha_{x})_{x\in
E},\vec{\beta}=(\beta_{x})_{x\in E}\in\mathbb{C}^{4}$. Since (233) and (234)
are equivalent to the statement that

$\displaystyle\|\vec{\alpha}\|^{2}+\|\vec{\beta}\|^{2}$ $\displaystyle=1$
(238) $\displaystyle\langle\vec{\alpha},\vec{\beta}\rangle$ $\displaystyle=0,$
(239)

we obtain

$\displaystyle\left\|\sqrt{f}\vec{\alpha}+\sqrt{1-f}\vec{\beta}\right\|^{2}$
$\displaystyle=f\|\vec{\alpha}\|^{2}+(1-f)\|\vec{\beta}\|^{2}$ (240)
$\displaystyle\leq\max\\{f,1-f\\}$ (241) $\displaystyle=f\ $ (242)

by the Pythagorean theorem in $\mathbb{C}^{4}$ (using (239)) and by maximizing
over $(\alpha,\beta)$ satisfying (238). Inserting (242) into (237) results in
the upper bound (230) on $\|\Pi_{E}\Psi\|^{2}$.

It remains to prove the claimed identities (233) and (234). We argue that
these are a consequence of the fact that $\Psi_{g}$ is normalized and
Gaussian, respectively.

###### Proof.

Observe that by definition of the antiunitary $\theta$, we have

$\displaystyle\beta_{x}$ $\displaystyle=\langle x,\theta\Psi_{g}\rangle$ (243)
$\displaystyle=(-1)^{\vartheta(\overline{x})}\overline{\langle\overline{x},\Psi_{g}\rangle}\qquad\textrm{
for all }\qquad x\in E\ .$ (244)

In particular, this implies that

$\displaystyle|\beta_{x}|^{2}=|\langle\overline{x},\Psi_{g}\rangle|^{2}\qquad\textrm{
for every }\qquad x\in E\ .$ (245)

Eq. (233) now follows from the fact that $\Psi_{g}$ is normalized and
positive-parity: we have

$\displaystyle\sum_{x\in E}\left(|\alpha_{x}|^{2}+|\beta_{x}|^{2}\right)$
$\displaystyle=\sum_{x\in E}\left(|\langle
x,\Psi_{g}\rangle|^{2}+|\langle\overline{x},\Psi_{g}\rangle|^{2}\right)$ (246)
$\displaystyle=\|\Psi_{g}\|^{2}$ (247) $\displaystyle=1$ (248)

where we used the definition of $\alpha_{x}$ and (245) in the first step, and
the assumption $E\cup\overline{E}=\\{0,1\\}^{4}_{+}$ in the second identity.

Similarly, Eq. (234) is a consequence of the fact that $\Psi_{g}$ is Gaussian:
we have

$\displaystyle\sum_{x\in\\{0,1\\}_{+}^{4}}\overline{\alpha_{x}}\beta_{x}$
$\displaystyle=\sum_{x\in\\{0,1\\}_{+}^{4}}\langle\Psi_{g},x\rangle\langle
x,\theta\Psi_{g}\rangle$ (249)
$\displaystyle=\langle\Psi_{g},\theta\Psi_{g}\rangle$ (250) $\displaystyle=0$
(251)

where we used the definition of $\alpha_{x}$ and $\beta_{x}$ in the first
step, the fact that $\Psi_{g}\in\mathcal{H}_{+}^{4}$ in the second step, and
the characterization of Gaussianity from Lemma 5.1 in the last identity. ∎

∎

Lemma 5.5 immediately implies the following expression for the fermionic
Gaussian fidelity. We note that a more general expression for the “Gaussian
fidelity” of a mixed state has previously been obtained in [58]. The proof for
pure states given here is more elementary and illustrates the use of Lemma
5.5.

######  Theorem 5.6 (Fermionic Gaussian fidelity for $4$-mode pure states [58,
56]).

Let $\Psi\in\mathcal{H}_{+}^{4}$ be a unit vector. Let $f(\Psi)\in[1/2,1]$ be
defined as in Theorem 5.4. Then

$\displaystyle F_{\mathcal{G}^{+}_{4}}(\Psi)$ $\displaystyle=f(\Psi)\ .$ (252)

###### Proof.

Let $f=f(\Psi)$, $\varphi\in[0,2\pi)$ and $\Psi_{g}\in\mathcal{G}^{+}_{4}$ be
as in Theorem 5.4. Then we have

$\displaystyle F_{\mathcal{G}^{+}_{4}}(\Psi)$
$\displaystyle\geq\left|\langle\Psi_{g},e^{i\varphi}\left(\sqrt{f}\Psi_{g}+f\theta\Psi_{g}\right)\rangle\right|^{2}$
(253) $\displaystyle=f\ $ (254)

since $\theta\Psi_{g}$ is orthogonal to $\Psi_{g}$. It thus suffices to show
the upper bound

$\displaystyle F_{\mathcal{G}^{+}_{4}}(\Psi)$ $\displaystyle\leq f\ .$ (255)

Let $\Phi_{g}\in\mathcal{G}^{+}_{4}$ be an arbitrary positive-parity Gaussian
pure state. Then there is a Gaussian unitary $U=U_{R}$ with $R\in SO(2n)$ and
a phase $\mu\in[0,2\pi)$ such that $\Phi_{g}=e^{i\mu}U\ket{0_{F}}$. We will
use any subset $E\subset\\{0,1\\}_{+}^{4}$ of even-weight strings as in Lemma
5.5 with the additional property that $0000\in E$, e.g.,
$E=\\{0000,1100,1010,1001\\}$. Then $\ket{0_{F}}=\Pi_{E}\ket{0_{F}}$ is in the
image of $\Pi_{E}$. It follows that

$\displaystyle\left|\langle\Phi_{g},\Psi\rangle\right|$
$\displaystyle=\left|\langle 0_{F},U^{\dagger}\Psi\rangle\right|$ (256)
$\displaystyle=\left|\langle\Pi_{E}0_{F},U^{\dagger}\Psi\rangle\right|$ (257)
$\displaystyle=\left|\langle 0_{F},\Pi_{E}U^{\dagger}\Psi\rangle\right|$ (258)
$\displaystyle\leq\left\|\Pi_{E}U^{\dagger}\Psi\right\|$ (259)
$\displaystyle\leq\sqrt{f(U^{\dagger}\Psi)}\ ,$ (260)

where we used the Cauchy-Schwarz inequality in the penultimate step, and Lemma
5.5 applied to the state $U^{\dagger}\Psi$. Since
$\Phi_{g}\in\mathcal{G}^{+}_{4}$ was arbitrary, the claimed inequality (255)
follows by taking the square and using that $f(U^{\dagger}\Psi)=f(\Psi)$, see
Eq. (221) of Theorem 5.4. ∎

Combining Lemma 5.5 with Theorem 5.6 yields the following statement, which
directly relates the weight of a state on certain subspaces to the fermionic
Gaussian fidelity. It will be our main technical tool in the following.

######  Corollary 5.7.

Let $E\subset\\{0,1\\}_{+}^{4}$, $|E|=4$ be a subset of even-weight strings
such that $E\cup\overline{E}=\\{0,1\\}^{4}_{+}$. Then

$\displaystyle\left\|\Pi_{E}\Psi\right\|^{2}\leq
F_{\mathcal{G}^{+}_{4}}(\Psi)\qquad\textrm{ for any unit vector
}\qquad\Psi\in\mathcal{H}_{+}^{4}\ .$ (261)

###  5.3 Multiplicativity of the Gaussian fidelity for 4-fermion states

Here we prove that the fermionic Gaussian fidelity is multiplicative for
$4$-fermion states (see Theorem 5.10). We use two intermediate results stated
as Lemmas 5.8 and 5.9. In Lemma 5.8, we bound the overlap of a tensor product
of two (arbitrary) positive-parity pure states $\Psi_{A},\Psi_{B}$ with a
state $\Phi$ written as a Schmidt decomposition of a bipartite fermionic pure
state. In Theorem 5.10, this result is used to bound the fermionic Gaussian
fidelity. More specifically, Lemma 5.9 is used to upper bound the Schmidt
coefficients, giving the multiplicativity result for the Gaussian fidelity.

######  Lemma 5.8.

Let $\\{m_{x}\\}_{x\in\\{0,1\\}^{4}}\subset\mathbb{C}$ be arbitrary. Define

$\displaystyle\ket{\Phi}=\sum_{x\in\\{0,1\\}^{4}}m_{x}\ket{x,x}\in\mathcal{H}^{8}_{+}\
.$ (262)

Let $E\subset\\{0,1\\}_{+}^{4}$, $|E|=4$ be a subset of even-weight strings
such that $E\cup\overline{E}=\\{0,1\\}^{4}_{+}$. Then

$\displaystyle|\langle\Phi,\Psi_{A}\tilde{\otimes}\Psi_{B}\rangle|^{2}$
$\displaystyle\leq
F_{\mathcal{G}^{+}_{4}}(\Psi_{A})F_{\mathcal{G}^{+}_{4}}(\Psi_{B})\left(\max_{x\in
E}|m_{x}|+\max_{y\in\overline{E}}|m_{y}|\right)^{2}$ (263)

for all states $\Psi_{A},\Psi_{B}\in\mathcal{H}_{+}^{4}$.

###### Proof.

Because $\Psi_{A}$ and $\Psi_{B}$ are supported on $\mathcal{H}_{+}^{4}$ by
assumption, we have by Eq. (74)

$\displaystyle\langle\Phi,\Psi_{A}\tilde{\otimes}\Psi_{B}\rangle$
$\displaystyle=\sum_{x\in\\{0,1\\}_{+}^{4}}m_{x}(-1)^{|x|}\langle
x,\Psi_{A}\rangle\langle x,\Psi_{B}\rangle$ (264)
$\displaystyle=\sum_{x\in\\{0,1\\}_{+}^{4}}m_{x}e^{i\nu_{x}}\langle\Psi_{A},x\rangle\langle
x,\Psi_{B}\rangle$ (265)

where $\nu_{x}$ is defined by the identity

$\displaystyle(-1)^{|x|}\langle x,\Psi_{A}\rangle$
$\displaystyle=e^{i\nu_{x}}\langle\Psi_{A},x\rangle\qquad\textrm{ for }\qquad
x\in\\{0,1\\}_{+}^{4}\ .$ (266)

Defining the operator

$\displaystyle M_{\Omega}=\sum_{x\in\Omega}m_{x}e^{i\nu_{x}}|x\rangle\langle
x|$ (267)

for any subset $\Omega\subset\\{0,1\\}_{+}^{4}$, it follows from Eq. (265)
that

$\displaystyle\langle\Phi,\Psi_{A}\tilde{\otimes}\Psi_{B}\rangle$
$\displaystyle=\langle\Psi_{A},M_{E}\Psi_{B}\rangle+\langle\Psi_{A},M_{\overline{E}}\Psi_{B}\rangle\
.$ (268)

Since $M_{E}$ is supported on $\mathsf{span}\\{\ket{x}\\}_{x\in E}$, we have
$\langle\Psi_{A},M_{E}\Psi_{B}\rangle=\langle\Pi_{E}\Psi_{A},M_{E}\Pi_{E}\Psi_{B}\rangle$.
With the Cauchy-Schwarz inequality and the definition of the operator norm
$\|M_{E}\|$ we thus get

$\displaystyle\left|\langle\Psi_{A},M_{E}\Psi_{B}\rangle\right|$
$\displaystyle\leq\|\Pi_{E}\Psi_{A}\|\cdot\|\Psi_{E}\Psi_{B}\|\cdot\|M_{E}\|$
(269)
$\displaystyle\leq\sqrt{F_{\mathcal{G}^{+}_{4}}(\Psi_{A})}\cdot\sqrt{F_{\mathcal{G}^{+}_{4}}(\Psi_{B})}\cdot\|M_{E}\|\
,$ (270)

where we applied Corollary 5.7. Identical reasoning applies to $\overline{E}$
and yields the inequality

$\displaystyle\left|\langle\Psi_{A},M_{\overline{E}}\Psi_{B}\rangle\right|$
$\displaystyle\leq\sqrt{F_{\mathcal{G}^{+}_{4}}(\Psi_{A})}\cdot\sqrt{F_{\mathcal{G}^{+}_{4}}(\Psi_{B})}\cdot\|M_{\overline{E}}\|\
.$ (271)

Combining Eqs. (270), (271) with Eq. (268), we conclude that

$\displaystyle\left|\langle\Phi,\Psi_{A}\tilde{\otimes}\Psi_{B}\rangle\right|$
$\displaystyle\leq\left|\langle\Psi_{A},M_{E}\Psi_{B}\rangle\right|+\left|\langle\Psi_{A},M_{\overline{E}}\Psi_{B}\rangle\right|$
(272)
$\displaystyle\leq\sqrt{F_{\mathcal{G}^{+}_{4}}(\Psi_{A})F_{\mathcal{G}^{+}_{4}}(\Psi_{B})}\left(\|M_{E}\|+\|M_{\overline{E}}\|\right)\
.$ (273)

Taking the square and observing that

$\displaystyle\|M_{\Omega}\|$
$\displaystyle=\max_{x\in\Omega}|m_{x}e^{i\nu_{x}}|=\max_{x\in\Omega}|m_{x}|\
$ (274)

gives the claim. ∎

The following lemma will be useful to prove the main theorem.

######  Lemma 5.9.

The function

$\displaystyle f(\theta,x)$
$\displaystyle=\prod_{j=1}^{4}(\cos\theta_{j})^{1-x_{j}}(\sin\theta_{j})^{x_{j}}\qquad\textrm{
for }\qquad\theta=(\theta_{1},\ldots,\theta_{4})\in\mathbb{R}^{4}\textrm{ and
}x\in\\{0,1\\}^{4}$ (275)

satisfies

$\displaystyle|f(\theta,x)|+|f(\theta,y)|$ $\displaystyle\leq 1$ (276)

for all $\theta\in\mathbb{R}^{4}$ and $x,y\in\\{0,1\\}^{4}$ with $x,y$ even-
weight and $x\neq y$.

###### Proof.

Because $x,y$ have even and different weight, it suffices to consider two
cases, namely with $|x-y|\in\\{2,4\\}$.

Consider first the case where $|x-y|=2$. Without loss of generality, assume
that $(x_{1},x_{2})=(y_{1},y_{2})$, $x_{3}\neq y_{3}$, $x_{4}\neq y_{4}$.
Since translating $\theta$ by $-\pi/2$ interchanges $\left|\sin\theta\right|$
and $\left|\cos\theta\right|$, it suffices to show the claim for $x=(0,0,0,0)$
and $y=(0,0,1,1)$. In this case we have

$\displaystyle|f(\theta,x)|+|f(\theta,y)|$
$\displaystyle=\left|\cos\theta_{1}\cos\theta_{2}\cos\theta_{3}\cos\theta_{4}\right|+\left|\cos\theta_{1}\cos\theta_{2}\sin\theta_{3}\sin\theta_{4}\right|$
(277)
$\displaystyle=\left|\cos\theta_{1}\cos\theta_{2}\right|\cdot\left(|\cos\theta_{3}\cos\theta_{4}|+|\sin\theta_{3}\sin\theta_{4}|\right)$
(278) $\displaystyle\leq\left|\cos\theta_{1}\cos\theta_{2}\right|$ (279)
$\displaystyle\leq 1\ ,$ (280)

where the first inequality follows from the Cauchy-Schwarz inequality in
$\mathbb{R}^{2}$. Since $\theta\in\mathbb{R}^{4}$ was arbitrary, this
concludes the proof for $|x-y|=2$.

The proof for $|x-y|=4$, i.e., $y=\overline{x}$ proceeds similarly. Again it
suffices to show the claim for $x=(0,0,0,0)$. In this case

$\displaystyle|f(\theta,x)|+|f(\theta,y)|$
$\displaystyle=|\cos\theta_{1}\cos\theta_{2}\cos\theta_{3}\cos\theta_{4}|+|\sin\theta_{1}\sin\theta_{2}\sin\theta_{3}\sin\theta_{4}|$
(281)
$\displaystyle\leq|\cos\theta_{1}\cos\theta_{2}|+|\sin\theta_{1}\sin\theta_{2}|$
(282) $\displaystyle\leq 1\ ,$ (283)

where the first inequality follows from $|\cos\theta_{3}\cos\theta_{4}|\leq
1~{}$ and $|\sin\theta_{3}\sin\theta_{4}|\leq 1~{}$ and the second one from
the Cauchy-Schwarz inequality. ∎

######  Theorem 5.10 (Multiplicativity of the fermionic Gaussian fidelity for
4-mode pure states.).

Let $\mathcal{H}^{n}_{+}$ be the set of pure $n$-fermion states with positive
parity and let $\mathcal{G}_{n}^{+}$ be the set of pure $n$-fermion Gaussian
states with positive parity. We have that

$\displaystyle
F_{\mathcal{G}_{8}}(\Psi_{A}\tilde{\otimes}\Psi_{B})=F_{\mathcal{G}_{4}}(\Psi_{A})F_{\mathcal{G}_{4}}(\Psi_{B})\qquad\textrm{
for all }\qquad\Psi_{A},\Psi_{B}\in\mathcal{H}^{4}_{+}\ .$ (284)

###### Proof.

We first observe that $\mathcal{H}^{8}$ is a direct sum of the four spaces
$\mathcal{H}^{4}_{+}\tilde{\otimes}\mathcal{H}^{4}_{+}$,
$\mathcal{H}^{4}_{+}\tilde{\otimes}\mathcal{H}^{4}_{-}$,
$\mathcal{H}^{4}_{-}\tilde{\otimes}\mathcal{H}^{4}_{+}$ and
$\mathcal{H}^{4}_{-}\tilde{\otimes}\mathcal{H}^{4}_{-}$. This is because
states in these subspaces have different eigenvalues with respect to the
corresponding parity operators on the factors (interpreted as Majorana
monomials on $\mathcal{H}^{8}$ these are the monomials $c(1^{8}0^{8})$ and
$c(0^{8}1^{8})$). It follows immediately that the overlap with a state of the
form
$\Psi_{A}\tilde{\otimes}\Psi_{B}\in\mathcal{H}^{4}_{+}\tilde{\otimes}\mathcal{H}^{4}_{+}$
is maximized for a decomposition into states belonging to
$\mathcal{H}^{4}_{+}\tilde{\otimes}\mathcal{H}^{4}_{+}$ only. In particular,
it follows that

$\displaystyle F_{\mathcal{G}_{8}}(\Psi_{A}\tilde{\otimes}\Psi_{B})$
$\displaystyle=F_{\mathcal{G}_{8}^{+}}(\Psi_{A}\tilde{\otimes}\Psi_{B})\qquad\textrm{
for all }\qquad\Psi_{A},\Psi_{B}\in\mathcal{H}^{4}_{+}\ $ (285)

and by the same reasoning, we have

$\displaystyle
F_{\mathcal{G}_{4}}(\Psi)=F_{\mathcal{G}_{4}^{+}}(\Psi)\qquad\textrm{ for all
}\qquad\Psi\in\mathcal{H}^{4}_{+}\ .$ (286)

We conclude that it suffices to show that

$\displaystyle
F_{\mathcal{G}_{8}^{+}}(\Psi_{A}\tilde{\otimes}\Psi_{B})=F_{\mathcal{G}_{4}^{+}}(\Psi_{A})F_{\mathcal{G}_{4}^{+}}(\Psi_{B})\qquad\textrm{
for all }\qquad\Psi_{A},\Psi_{B}\in\mathcal{H}^{4}_{+}\ .$ (287)

Let $\Psi_{A},\Psi_{B}\in\mathcal{H}^{4}_{+}$ be arbitrary. The inequality
$F_{\mathcal{G}_{8}^{+}}(\Psi_{A}\tilde{\otimes}\Psi_{B})\geq
F_{\mathcal{G}_{4}^{+}}(\Psi_{A})F_{\mathcal{G}_{4}^{+}}(\Psi_{B})$ follows
trivially from the definition of fermionic Gaussian fidelity in Eq. (4),
because
$\mathcal{G}_{4}^{+}\tilde{\otimes}\mathcal{G}_{4}^{+}\subseteq\mathcal{G}_{8}^{+}$.

The inequality $F_{\mathcal{G}_{8}^{+}}(\Psi_{A}\tilde{\otimes}\Psi_{B})\leq
F_{\mathcal{G}_{4}^{+}}(\Psi_{A})F_{\mathcal{G}_{4}^{+}}(\Psi_{B})$ is a
consequence of the Schmidt decomposition for fermionic states put forward in
Ref. [50] and of Lemmas 5.8 and 5.9. According to Ref. [50], an arbitrary pure
fermionic state $\Phi\in\mathcal{H}^{n}$ admits a Schmidt decomposition of the
form

$\displaystyle|\Phi\rangle=\sum_{x\in\\{0,1\\}^{n}}m_{x}|x,x\rangle$ (288)

with

$\displaystyle
m_{x}=\prod_{j=1}^{n}(\cos\theta_{j})^{1-x_{j}}(-\sin\theta_{j})^{x_{j}}\qquad\mathrm{with}\qquad\theta_{j}\in\mathbb{R}\text{
for }j\in[n]\ .$ (289)

With this definition of $m_{x}$ for $n=4$, an arbitrary state
$\Phi\in\mathcal{G}_{8}^{+}$ can be written as in Eq. (262) and the conditions
for Lemma 5.8 apply. We have

$\displaystyle F_{\mathcal{G}_{8}^{+}}(\Psi_{A}\tilde{\otimes}\Psi_{B})$
$\displaystyle=\max_{\Phi\in\mathcal{G}_{8}^{+}}|\langle\Phi,\Psi_{A}\otimes\Psi_{B}\rangle|$
(290) $\displaystyle\leq
F_{\mathcal{G}^{+}_{4}}(\Psi_{A})F_{\mathcal{G}^{+}_{4}}(\Psi_{B})\left(\max_{x\in
E}|m_{x}|+\max_{y\in\overline{E}}|m_{y}|\right)^{2}$ (291)

where $E\subset\\{0,1\\}^{4}_{+}$ with $|E|=4$ a subset of even weight strings
such that $E\cup\overline{E}=\\{0,1\\}^{4}_{+}$. Notice that $x,y$ have even-
weight and that $x\neq y$ because $E$ and $\overline{E}$ are disjoint sets.
Identifying $m_{x}$ (whose dependence on $\theta\in\mathbb{R}^{4}$ is
implicit) with $f(\theta,x)$ in Lemma 5.9 (apart from a minus sign that is not
relevant because we take the absolute value) we have

$\displaystyle F_{\mathcal{G}_{8}^{+}}(\Psi_{A}\tilde{\otimes}\Psi_{B})$
$\displaystyle\leq
F_{\mathcal{G}^{+}_{4}}(\Psi_{A})F_{\mathcal{G}^{+}_{4}}(\Psi_{B})\left(\max_{x\in
E}|m_{x}|+\max_{y\in\overline{E}}|m_{y}|\right)^{2}$ (292) $\displaystyle\leq
F_{\mathcal{G}^{+}_{4}}(\Psi_{A})F_{\mathcal{G}^{+}_{4}}(\Psi_{B})\ ,$ (293)

giving the claim. ∎

##  6 Multiplicativity of $\mathcal{D}$-fidelity implies that of
$\mathcal{D}$-extent

In this section, we show that multiplicativity of the $\mathcal{D}$-fidelity
implies multiplicativity of the $\mathcal{D}$-extent. In Section 6.1 we prove
this for finite dictionaries: This follows immediately from the fact that
$F_{\mathcal{D}}(\Psi)$ and $\xi_{\mathcal{D}}(\Psi)$ are related by (convex
programming) duality. In Section 6.2, we extend this results for infinite,
i.e., continuously parameterized dictionaries. We achieve this extension by
using (finite) $\epsilon$-nets for the set of Gaussian states. Similar
approaches have been applied in the signal processing context, see e.g., the
work [63], which shows how to approximately solve atomic norm minimization
problems for sparse recovery when the parameters indexing the dictionary lie
in a small-dimensional space.

###  6.1 Multiplicativity for finite dictionaries

We will restrict our attention to finite dictionaries in this section. For
$|\mathcal{D}|<\infty$, the $\mathcal{D}$-fidelity is related to the dual
formulation of the $\mathcal{D}$-extent as (see [19, Eq. (3.2)] and [18,
Theorem 4])

$\displaystyle\xi_{\mathcal{D}}(\Psi)$
$\displaystyle=\max_{y\in\mathcal{H}:F_{\mathcal{D}}(y)\leq
1}|\langle\Psi,y\rangle|^{2}\ \ .$ (294)

Let $\mathcal{H}_{1},\mathcal{H}_{2}$ and $\mathcal{H}_{3}$ be a triple of
Hilbert spaces. Let $\\{\mathcal{D}_{j}\\}_{j\in[3]}$ be a family of
dictionaries, where $\mathcal{D}_{j}\subset\mathcal{H}^{j}$ for $j\in[3]$. We
assume that

$\displaystyle\mathcal{D}_{1}\otimes\mathcal{D}_{2}\subseteq\mathcal{D}_{3}\
.$ (295)

We are interested in the following two properties:

$\displaystyle\mathsf{Mult}^{\xi}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$
$\displaystyle:\quad\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})=\xi_{\mathcal{D}_{1}}(\Psi_{1})\xi_{\mathcal{D}_{2}}(\Psi_{2})\quad\textrm{
for all }\quad\Psi_{j}\in\mathcal{H}_{j}\textrm{ for }j\in[2]\ $ (296)
$\displaystyle\mathsf{Mult}^{F}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$
$\displaystyle:\quad
F_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})=F_{\mathcal{D}_{1}}(\Psi_{1})F_{\mathcal{D}_{2}}(\Psi_{2})\quad\textrm{
for all }\quad\Psi_{j}\in\mathcal{H}_{j}\textrm{ for }j\in[2]\ .$ (297)

As an important example, let $n_{j}\in\mathbb{N}$ for $j\in[2]$,
$n_{3}=n_{1}+n_{2}$, $\mathcal{H}_{j}=(\mathbb{C}^{2})^{\otimes n_{j}}$ and
let $\mathsf{STAB}_{n}$ be the set of stabilizer states on
$(\mathbb{C}^{2})^{\otimes n}$. Then
$\mathsf{Mult}^{\xi}(\mathsf{STAB}_{n_{1}},\mathsf{STAB}_{n_{1}},\mathsf{STAB}_{n_{1}+n_{2}})$
does not hold for certain (large) choices of $n_{1}$ and $n_{2}$ [19]. On the
other hand, for $n_{1},n_{2}\leq 3$ the multiplicativity property
$\mathsf{Mult}^{\xi}(\mathsf{STAB}_{n_{1}},\mathsf{STAB}_{n_{1}},\mathsf{STAB}_{n_{1}+n_{2}})$
holds (see Ref. [18, Proposition 1]). This was shown using that the stabilizer
fidelity is multiplicative, i.e.,
$\mathsf{Mult}^{F}(\mathsf{STAB}_{n_{1}},\mathsf{STAB}_{n_{1}},\mathsf{STAB}_{n_{1}+n_{2}})$.

We claim the property
$\mathsf{Mult}^{F}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$ implies
property
$\mathsf{Mult}^{\xi}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$.

######  Theorem 6.1.

Property $\mathsf{Mult}^{F}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$
implies property
$\mathsf{Mult}^{\xi}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$.

###### Proof.

We clearly have

$\displaystyle\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})$
$\displaystyle\leq\xi_{\mathcal{D}_{1}}(\Psi_{1})\xi_{\mathcal{D}_{2}}(\Psi_{2})$
(298)

for all $\Psi_{1}\in\mathcal{H}_{1}$ and $\Psi_{2}\in\mathcal{H}_{2}$ because
of property (295) of the dictionaries $\\{\mathcal{D}_{j}\\}_{j=1}^{3}$ and
the definition (8) of $\xi_{\mathcal{D}}$. To show the converse inequality,
assume that $y_{1}\in\mathcal{H}_{1}$, $y_{2}\in\mathcal{H}_{2}$ are such that

$\displaystyle F_{\mathcal{D}_{j}}(y_{j})\leq 1\quad\textrm{ and
}\quad\xi_{\mathcal{D}_{j}}(\Psi_{j})=|\langle\Psi_{j},y_{j}\rangle|^{2}\quad\textrm{
for }\quad j\in[2]\ .$ (299)

Then

$\displaystyle F_{\mathcal{D}_{3}}(y_{1}\otimes y_{2})$
$\displaystyle=F_{\mathcal{D}_{1}}(y_{1})F_{\mathcal{D}_{2}}(y_{2})\leq 1$
(300)

where we used the assumption that property
$\mathsf{Mult}^{F}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$ holds to
obtain the equality. This implies that $y_{1}\otimes y_{2}$ is a feasible
point of the dual program for the quantity
$\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})$, see Eq. (294). Thus

$\displaystyle\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})$
$\displaystyle\geq|\langle\Psi_{1}\otimes\Psi_{2},y_{1}\otimes
y_{2}\rangle|^{2}$ (301)
$\displaystyle=|\langle\Psi_{1},y_{1}\rangle|^{2}\cdot|\langle\Psi_{2},y_{2}\rangle|^{2}$
(302)
$\displaystyle=\xi_{\mathcal{D}_{1}}(\Psi_{1})\xi_{\mathcal{D}_{2}}(\Psi_{2})\
.$ (303)

Expression (303) together with Eq. (298) gives the claim. ∎

###  6.2 Multiplicativity for infinite dictionaries

In this section, we extend the results of Section 6.2 to dictionaries
$\mathcal{D}$ that may contain infinitely many elements. Our strategy is to
use an $\epsilon$-net for $\mathcal{D}\in\mathcal{H}$ with a finite number of
elements, we denote by $\mathcal{D}^{\epsilon}$. We relate the extent and
fidelity with respect to the dictionary $\mathcal{D}$ to the extent and
fidelity with respect to its net $\mathcal{D}^{\epsilon}$ (see Lemmas 6.2 and
6.3) to prove that multiplicativity of the $\mathcal{D}$-fidelity implies
multiplicativity of the $\mathcal{D}$-extent in Theorem 6.6. This result is a
generalization of Theorem 6.1 (that considered finite dictionaries) for
(possibly) infinite dictionaries.

We will make use of the notion of $\epsilon$-net to replace our infinite set
$\mathcal{D}$ by a finite set $\mathcal{D}^{\epsilon}$. Let
$\|\Psi\|=\sqrt{\langle\Psi,\Psi\rangle}$ for $\Psi\in\mathcal{H}$ denote the
norm on $\mathcal{H}$. Let $\mathcal{D}\subset\mathcal{H}$ and let
$\epsilon>0$. Then a set $\mathcal{D}^{\epsilon}\subset\mathcal{H}$ is called
an $\epsilon$-net for $\mathcal{D}$ if for any $\Psi\in\mathcal{D}$ there is
some $\Phi\in\mathcal{D}^{\epsilon}$ such that $\|\Phi-\Psi\|\leq\epsilon$.

We are interested in the case where for every $\epsilon>0$ there is a finite
$\epsilon$-net $\mathcal{D}^{\epsilon}$ for $\mathcal{D}$, with the additional
property that $\mathcal{D}^{\epsilon}\subset\mathcal{D}$, i.e., the net
consists of elements of $\mathcal{D}$. A sufficient condition for this being
the case is that the subset $\mathcal{D}\subset\mathcal{H}$ is compact.

######  Lemma 6.2.

Let $\mathcal{D}\subset\mathcal{H}$ be a set of states. Assume that there is a
finite $\epsilon$-net $\mathcal{D}^{\epsilon}$ for $\mathcal{D}$ such that
$\mathcal{D}^{\epsilon}\subset\mathcal{D}$, for some $\epsilon>0$. Assume
further that $\mathcal{D}^{\epsilon}$ contains an orthonormal basis of
$\mathcal{H}$. Let $d$ be the dimension of $\mathcal{H}$. Then

$\displaystyle\xi_{\mathcal{D}}(\Psi)$
$\displaystyle\leq\xi_{\mathcal{D}^{\epsilon}}(\Psi)\leq\xi_{\mathcal{D}}(\Psi)\left(1+\sqrt{d}\epsilon\right)^{2}\qquad\textrm{
for all }\qquad\Psi\in\mathcal{H}\ .$ (304)

###### Proof.

The first inequality follows immediately from the definition of
$\xi_{\mathcal{D}}$ and from the assumption that
$\mathcal{D}^{\epsilon}\subset\mathcal{D}$.

To prove the second inequality, let $\Psi\in\mathcal{H}$ be arbitrary. By
definition of $\xi_{\mathcal{D}}(\Psi)$ as an infimum, we have the following:
For every $m\in\mathbb{N}$, there exist $N(m)\in\mathbb{N}$,
$\\{\varphi_{j}(m)\\}_{j=1}^{N(m)}\subset\mathcal{D}$ and
$\\{c_{j}(m)\\}_{j=1}^{N(m)}\subset\mathbb{C}$ such that

$\displaystyle\Psi=\sum_{j=1}^{N(m)}c_{j}(m)\varphi_{j}(m)$ (305)

and

$\displaystyle\|c(m)\|_{1}$
$\displaystyle<\sqrt{\xi_{\mathcal{D}}(\Psi)}+\frac{1}{m}\ .$ (306)

Furthermore, we have

$\displaystyle\|c(m)\|_{1}^{2}$
$\displaystyle\geq\xi_{\mathcal{D}}(\Psi)\quad\textrm{ for any
}\quad{N(m)}\in\mathbb{N},\\{s_{j}\\}_{j=1}^{N(m)}\subset\mathcal{D}\quad\textrm{
and }\quad\\{c_{j}(m)\\}_{j=1}^{N(m)}\subset\mathbb{C}\ .$ (307)

Fix such an $m\in\mathbb{N}$. Since $\mathcal{D}^{\epsilon}$ is an
$\epsilon$-net for $\mathcal{D}$, there is, for every $j\in[N(m)]$, an element
$\varphi_{j}^{\epsilon}(m)\in\mathcal{D}^{\epsilon}$ and
$\delta_{j}(m)\in\mathcal{H}$ such that

$\displaystyle\varphi_{j}(m)$
$\displaystyle=\varphi_{j}^{\epsilon}(m)+\delta_{j}(m)\qquad\textrm{ and
}\qquad\|\delta_{j}(m)\|\leq\epsilon\ .$ (308)

It follows that

$\displaystyle\Psi$
$\displaystyle=\left(\sum_{j=1}^{N(m)}c_{j}(m)\varphi_{j}^{\epsilon}(m)\right)+\delta(m)\qquad\textrm{
where }\qquad\delta(m)=\sum_{j=1}^{N(m)}c_{j}(m)\delta_{j}(m)\ .$ (309)

By the triangle inequality we have

$\displaystyle\|\delta(m)\|$
$\displaystyle\leq\epsilon\sum_{j=1}^{N(m)}|c_{j}(m)|=\epsilon\|c(m)\|_{1}\ .$
(310)

Suppose $\\{e_{k}\\}_{k=1}^{d}$ is an orthonormal basis contained in
$\mathcal{D}^{\epsilon}$. Then we can expand

$\displaystyle\delta(m)$ $\displaystyle=\sum_{k=1}^{d}\alpha_{k}(m)e_{k}\ ,$
(311)

and it follows from (310) and the Cauchy-Schwarz inequality that

$\displaystyle\|\alpha(m)\|_{1}$ $\displaystyle\leq\sqrt{d}\|\alpha(m)\|_{2}$
(312) $\displaystyle=\sqrt{d}\|\delta(m)\|_{2}$ by the Cauchy-Schwarz
inequality in $\mathbb{C}^{d}$ (313)
$\displaystyle\leq\sqrt{d}\epsilon\|c(m)\|_{1}\ $
$\displaystyle\qquad\text{because of Eq.~{}\eqref{eq:deltanormproperty}}\ .$
(314)

where $\|\cdot\|_{2}$ is the Euclidean norm in $\mathbb{C}^{d}$. In summary,
we have

$\displaystyle\Psi$
$\displaystyle=\sum_{j=1}^{N(m)}c_{j}(m)\varphi_{j}^{\epsilon}(m)+\sum_{k=1}^{d}\alpha_{k}(m)e_{k}\
,$ (315)

with $(c,\alpha)\in\mathbb{C}^{N(m)+d}$ satisfying

$\displaystyle\|(c,\alpha)\|_{1}$
$\displaystyle=\sum_{j=1}^{N(m)}|c_{j}(m)|+\sum_{k=1}^{d}|\alpha_{k}(m)|$
(316) $\displaystyle\leq\|c(m)\|_{1}\cdot\left(1+\sqrt{d}\epsilon\right)$
(317)
$\displaystyle\leq\left(\sqrt{\xi_{\mathcal{D}}(\Psi)}+\frac{1}{m}\right)\cdot\left(1+\sqrt{d}\epsilon\right)\
,$ (318)

by (306). Because $\varphi_{j}^{\epsilon}(m)\in\mathcal{D}^{\epsilon}$ for
every $j\in[N(m)]$ and $e_{k}\in\mathcal{D}^{\epsilon}$ for every $k\in[d]$,
this implies that

$\displaystyle\sqrt{\xi_{\mathcal{D}^{\epsilon}}(\Psi)}\leq\left(\sqrt{\xi_{\mathcal{D}}(\Psi)}+\frac{1}{m}\right)\cdot\left(1+\sqrt{d}\epsilon\right)\
.$ (319)

Since $m\in\mathbb{N}$ was arbitrary, we can take the limit
$m\rightarrow\infty$ and obtain

$\displaystyle\sqrt{\xi_{\mathcal{D}^{\epsilon}}(\Psi)}\leq\sqrt{\xi_{\mathcal{D}}(\Psi)}\cdot\left(1+\sqrt{d}\epsilon\right)\
.$ (320)

This gives the claim. ∎

######  Lemma 6.3.

Let $\mathcal{D}^{\epsilon}\subset\mathcal{D}\subset\mathcal{H}$ be as in
Lemma 6.2. Then we have

$\displaystyle\sqrt{F_{\mathcal{D}^{\epsilon}}(\Psi)}$
$\displaystyle\leq\sqrt{F_{\mathcal{D}}(\Psi)}\leq\sqrt{F_{\mathcal{D}^{\epsilon}}(\Psi)}+\|\Psi\|\cdot\epsilon\qquad\textrm{
for all }\qquad\Psi\in\mathcal{H}\ .$ (321)

###### Proof.

The first inequality follows trivially from the definitions using
$\mathcal{D}^{\epsilon}\subset\mathcal{D}$. For the second inequality, let
$\Psi\in\mathcal{H}$ be arbitrary. Let $\varphi\in\mathcal{D}$ be such that

$\displaystyle F_{\mathcal{D}}(\Psi)$
$\displaystyle=|\langle\varphi,\Psi\rangle|^{2}\ .$ (322)

Then (by the fact that $\mathcal{D}^{\epsilon}$ is an $\epsilon$-net for
$\mathcal{D}$ and $\varphi\in\mathcal{D}$) there is an element
$\varphi^{\epsilon}\in\mathcal{D}^{\epsilon}$ and $\delta\in\mathcal{H}$ such
that

$\displaystyle\varphi$ $\displaystyle=\varphi^{\epsilon}+\delta\qquad\textrm{
where }\qquad\|\delta\|\leq\epsilon\ .$ (323)

It follows that

$\displaystyle\sqrt{F_{\mathcal{D}}(\Psi)}$
$\displaystyle=|\langle\varphi,\Psi\rangle|$ (324)
$\displaystyle\leq|\langle\varphi^{\epsilon},\Psi\rangle|+|\langle\delta,\Psi\rangle|$
(325)
$\displaystyle\leq\sqrt{F_{\mathcal{D}^{\epsilon}}(\Psi)}+|\langle\delta,\Psi\rangle|$
(326)
$\displaystyle\leq\sqrt{F_{\mathcal{D}^{\epsilon}}(\Psi)}+\|\Psi\|\cdot\|\delta\|$
(327)
$\displaystyle\leq\sqrt{F_{\mathcal{D}^{\epsilon}}(\Psi)}+\|\Psi\|\cdot\epsilon$
(328)

where we used the definition of $F_{\mathcal{D}^{\epsilon}}(\Psi)$ and the
Cauchy-Schwarz inequality (in the penultimate step). The claim follows. ∎

######  Lemma 6.4.

We have

$\displaystyle\|y\|^{2}\leq d^{2}\cdot
F_{\mathcal{D}^{\epsilon}}(y)\qquad\textrm{ for every }\qquad y\in\mathcal{H}\
,$ (329)

where $d=\dim\mathcal{H}$.

###### Proof.

Let $\\{e_{k}\\}_{k=1}^{d}$ be an orthonormal basis contained in
$\mathcal{D}^{\epsilon}$. Then

$\displaystyle|\langle e_{k},y\rangle|^{2}$ $\displaystyle\leq
F_{\mathcal{D}^{\epsilon}}(y)\qquad\textrm{ for every }\qquad k\in[d]$ (330)

because $e_{k}\in\mathcal{D}^{\epsilon}$ for every $k\in[n]$. We have

$\displaystyle\|y\|$ $\displaystyle=\left(\sum_{k=1}^{d}|\langle
e_{k},y\rangle|^{2}\right)^{1/2}$ (331)
$\displaystyle\leq\sum_{k=1}^{d}|\langle e_{k},y\rangle|$ (332)
$\displaystyle\leq d\cdot\sqrt{F_{\mathcal{D}^{\epsilon}}(y)}$ (333)

where we used that $\|v\|_{2}\leq\|v\|_{1}$ for $v\in\mathbb{C}^{d}$. The
claim follows. ∎

######  Lemma 6.5.

Let $d_{j}=\dim\mathcal{H}_{j}$ for $j\in[2]$. Assuming that the property
$\mathsf{Mult}^{F}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$ holds, we
have

$\displaystyle\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})$
$\displaystyle\geq
g(d_{1},d_{2},d_{3},\epsilon)\xi_{\mathcal{D}_{1}}(\Psi_{1})\xi_{\mathcal{D}_{2}}(\Psi_{2})$
(334)

for a function $g$ which satisfies

$\displaystyle\lim_{\epsilon\rightarrow 0}g(d_{1},d_{2},d_{3},\epsilon)=1\ .$
(335)

###### Proof.

Suppose $y_{j}\in\mathcal{H}_{j}$ for $j\in[2]$ is such that

$\displaystyle F_{\mathcal{D}_{j}^{\epsilon}}(y_{j})$ $\displaystyle\leq 1$
(336)

and

$\displaystyle\xi_{\mathcal{D}_{j}^{\epsilon}}(\Psi_{j})$
$\displaystyle=|\langle y_{j},\Psi_{j}\rangle|^{2}\ $ (337)

for $j\in[2]$. Such a pair $(y_{1},y_{2})$ exists since
$\mathcal{D}^{\epsilon}_{1}$ and $\mathcal{D}^{\epsilon}_{2}$ are finite sets
and the dual definition definition of the extent in terms of a maximum applies
(see Eq. (294)). Equation (336) implies that

$\displaystyle\|y_{j}\|$ $\displaystyle\leq d_{j}\qquad\textrm{ for }\qquad
j\in[2]\ ,$ (338)

see Lemma 6.4.

We have

$\displaystyle\sqrt{F_{\mathcal{D}_{3}^{\epsilon}}(y_{1}\otimes y_{2})}$
$\displaystyle\leq\sqrt{F_{\mathcal{D}_{3}}(y_{1}\otimes
y_{2})}\qquad\qquad\qquad\qquad\qquad\qquad\,\qquad\textrm{ by
Lemma~{}\ref{lem:fapproximate}}$ (339)
$\displaystyle=\sqrt{F_{\mathcal{D}_{1}}(y_{1})}\cdot\sqrt{F_{\mathcal{D}_{2}}(y_{2})}\,\,\textrm{
by the assumption
that~{}$\mathsf{Mult}^{F}(\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3})$
holds}$ (340)
$\displaystyle\leq\left(\sqrt{F_{\mathcal{D}^{\epsilon}_{1}}(y_{1})}+\|y_{1}\|\cdot\epsilon\right)\left(\sqrt{F_{\mathcal{D}^{\epsilon}_{2}}(y_{2})}+\|y_{2}\|\cdot\epsilon\right)\quad\textrm{
by Lemma~{}\ref{lem:fapproximate}}$ (341)
$\displaystyle\leq\sqrt{F_{\mathcal{D}^{\epsilon}_{1}}(y_{1})}\cdot\sqrt{F_{\mathcal{D}^{\epsilon}_{2}}(y_{2})}\cdot\left(1+d_{1}\cdot\epsilon\right)\left(1+d_{2}\cdot\epsilon\right)\,\quad\textrm{
by Lemma~{}\ref{lem:boundednessyops}}$ (342)
$\displaystyle\leq(1+d_{1}\epsilon)(1+d_{2}\epsilon)\quad\qquad\qquad\qquad\qquad\qquad\qquad\,\,\textrm{by
Eq.~{}\eqref{eq:conditionyjupperbound}}\ .$ (343)

Defining

$\displaystyle\tilde{y}_{j}$
$\displaystyle=\left(1+d_{j}\epsilon\right)^{-1}y_{j}\qquad\textrm{ for
}\qquad j\in[2]$ (344)

it follows from Eq. (343) that

$\displaystyle
F_{\mathcal{D}_{3}^{\epsilon}}(\tilde{y}_{1}\otimes\tilde{y}_{2})$
$\displaystyle\leq 1\ .$ (345)

By the dual formulation (294) of the quantity
$\xi_{\mathcal{D}^{\epsilon}_{3}}(\Psi_{1}\otimes\Psi_{2})$ this implies that

$\displaystyle\xi_{\mathcal{D}^{\epsilon}_{3}}(\Psi_{1}\otimes\Psi_{2})$
$\displaystyle\geq|\langle\tilde{y}_{1}\otimes\tilde{y}_{2},\Psi_{1}\otimes\Psi_{2}\rangle|^{2}$
(346)
$\displaystyle=\left(\prod_{j=1}^{2}\left(1+d_{j}\epsilon\right)^{-2}|\langle
y_{j},\Psi_{j}\rangle|^{2}\right)$ (347)
$\displaystyle=\left(\prod_{j=1}^{2}\left(1+d_{j}\epsilon\right)^{-2}\right)\xi_{\mathcal{D}^{\epsilon}_{1}}(\Psi_{1})\xi_{\mathcal{D}^{\epsilon}_{2}}(\Psi_{2})$
because of Eq. (337) (348)
$\displaystyle\geq\left(\prod_{j=1}^{2}\left(1+d_{j}\epsilon\right)^{-2}\right)\xi_{\mathcal{D}_{1}}(\Psi_{1})\xi_{\mathcal{D}_{2}}(\Psi_{2})$
according to Lemma 6.2. (349)

Lemma 6.2 also implies that

$\displaystyle\left(1+\sqrt{d_{3}}\epsilon\right)^{2}\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})\geq\xi_{\mathcal{D}^{\epsilon}_{3}}(\Psi_{1}\otimes\Psi_{2})\
.$ (350)

Combining Eqs. (349) and (350) gives the claim, with

$\displaystyle
g(d_{1},d_{2},d_{3}\epsilon)=\left(1+\sqrt{d_{3}}\epsilon\right)^{-2}\left(\prod_{j=1}^{2}\left(1+d_{j}\epsilon\right)^{-2}\right)\
.$ (351)

∎

######  Theorem 6.6.

Let $\mathcal{H}_{1},\mathcal{H}_{2}$ and $\mathcal{H}_{3}$ be a triple of
Hilbert spaces and let $\\{\mathcal{D}_{j}\\}_{j\in[3]}$ be a family of
dictionaries, where $\mathcal{D}_{j}\subset\mathcal{H}_{j}$. Assume that
$\mathcal{D}_{j}$ contains an orthonormal basis of $\mathcal{H}_{j}$, for
$j\in[3]$, and that
$\mathcal{D}_{1}\otimes\mathcal{D}_{2}\subset\mathcal{D}_{3}$. Assume further
that for any $\epsilon>0$ there is an $\epsilon$-net
$\mathcal{D}_{j}^{\epsilon}$ for $\mathcal{D}_{j}$ such that
$\mathcal{D}_{j}^{\epsilon}\subset\mathcal{D}_{j}$, i.e., the net consists of
elements of $\mathcal{D}_{j}$. Finally, assume that

$\displaystyle F_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})$
$\displaystyle=F_{\mathcal{D}_{1}}(\Psi)F_{\mathcal{D}_{2}}(\Psi_{2})\qquad\textrm{
for all }\qquad\Psi_{1}\in\mathcal{H}_{1},\Psi_{2}\in\mathcal{H}_{2}\ .$ (352)

Then

$\displaystyle\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})$
$\displaystyle=\xi_{\mathcal{D}_{1}}(\Psi_{1})\xi_{\mathcal{D}_{2}}(\Psi_{2})\qquad\textrm{
for all }\qquad\Psi_{1}\in\mathcal{H}_{1}\textrm{ and
}\Psi_{2}\in\mathcal{H}_{2}\ .$ (353)

###### Proof.

By (if necessary) replacing $\mathcal{D}^{\epsilon}_{j}$ by
$\mathcal{D}^{\epsilon}_{j}\cup\\{e^{(j)}_{k}\\}_{k=1}^{d}$, for $j\in[3]$,
where $\\{e^{(j)}_{k}\\}_{k=1}^{d}$ is an orthonormal basis of
$\mathcal{H}_{j}$ with $e^{(j)}_{k}\in\mathcal{D}_{j}$, for $n\in[3]$, we have
that each $\mathcal{D}^{\epsilon}_{j}$ is finite and contains an orthonormal
basis of the respective space. The inequality

$\displaystyle\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})$
$\displaystyle\geq\xi_{\mathcal{D}_{1}}(\Psi_{1})\xi_{\mathcal{D}_{2}}(\Psi_{2})\qquad\textrm{
for all }\qquad\Psi_{1}\in\mathcal{H}_{1}\textrm{ and
}\Psi_{2}\in\mathcal{H}_{2}\ $ (354)

now follows immediately from Lemma 6.5 by taking the limit
$\epsilon\rightarrow 0$. The converse inequality is trivial because
$\mathcal{D}_{1}\otimes\mathcal{D}_{2}\subset\mathcal{D}_{3}$. ∎

##  7 Multiplicativity of the Gaussian extent for four fermions

In this section we prove that the Gaussian extent is multiplicative for the
tensor product of any two 4-fermion pure states with positive parity.

######  Theorem 7.1 (Multiplicativity of the Gaussian extent for 4-fermion
pure states.).

Let $\mathcal{H}^{4}_{+}$ be the set of pure $4$-fermion states with positive
parity and let $\mathcal{G}_{n}$ be the set of Gaussian states on $n$
fermions. Then

$\displaystyle\xi_{\mathcal{G}_{8}}(\Psi_{A}\tilde{\otimes}\Psi_{B})=\xi_{\mathcal{G}_{4}}(\Psi_{A})\xi_{\mathcal{G}_{4}}(\Psi_{B})\qquad\textrm{
for all }\qquad\Psi_{A},\Psi_{B}\in\mathcal{H}^{4}_{+}\ .$ (355)

###### Proof.

Since the metaplectic representation defines a surjective, continuous map

$\displaystyle\begin{matrix}f:&[0,2\pi]\times
SO(2n)&\rightarrow&\mathcal{G}_{n}\\\
&(\varphi,R)&\mapsto&e^{i\varphi}U_{R}\ket{0_{F}}\end{matrix}$ (356)

from the compact set $[0,2\pi]\times SO(2n)$ to $\mathcal{G}_{n}$, the set
$\mathcal{G}_{n}\subset\mathcal{H}^{n}$ is compact. We also observe that the
occupation number states

$\displaystyle\left\\{\ket{x}\ |\ x\in\\{0,1\\}^{n}\right\\}$ (357)

form an orthonormal basis contained in $\mathcal{G}_{n}$. By compactness, we
conclude that for any $\epsilon>0$, there is a finite $\epsilon$-net
$\mathcal{G}_{n}^{\epsilon}\subset\mathcal{G}_{n}$ consisting of Gaussian
states and containing an orthonormal basis of $\mathcal{H}^{n}$. Finally, we
note that we also have the inclusion
$\mathcal{G}_{n_{1}}\otimes\mathcal{G}_{n_{2}}\subset\mathcal{G}_{n_{1}+n_{2}}$
for $n_{1},n_{2}\in\mathbb{N}$ arbitrary.

Let us now specialize to $n_{1}=n_{2}=4$. In this case, we have
multiplicativity of the fermionic Gaussian fidelity by Theorem 5.10. In
particular, the conditions for Theorem 6.6 apply and the claim follows. ∎

## Acknowledgements

BD and RK gratefully acknowledge support by the European Research Council
under grant agreement no. 101001976 (project EQUIPTNT).

##  Appendix A Alternative Gaussianity condition for $4$-fermion states

In the following, we prove Lemma 5.1.

###### Proof.

Consider $\Psi\in\mathcal{H}_{+}^{4}$. We will show that
$\langle\Psi,\theta\Psi\rangle=0$ is equivalent to
$\Lambda\left(\ket{\Psi}\otimes\ket{\Psi}\right)=0$. According to Lemma 2.3,
this is a sufficient and necessary condition for $\Psi$ to be Gaussian.

Let $E\subset\\{0,1\\}_{+}^{4}$, $|E|=4$ be a subset of even-weight strings
such that $E\cup\overline{E}=\\{0,1\\}^{4}_{+}$. We first compute
$\langle\Psi,\theta\Psi\rangle$. We have

$\displaystyle\langle\Psi,\theta\Psi\rangle$
$\displaystyle=\sum_{x,y\in\\{0,1\\}^{4}_{+}}\langle\Psi,x\rangle\bra{x}\theta\left(\ket{y}\langle
y,\Psi\rangle\right)$ (358)
$\displaystyle=\sum_{x,y\in\\{0,1\\}^{4}_{+}}\langle\Psi,x\rangle\langle
x,\theta y\rangle\overline{\langle y,\Psi\rangle}$ (359)
$\displaystyle=\sum_{x,y\in\\{0,1\\}^{4}_{+}}(-1)^{\vartheta(y)}\langle
x,\overline{y}\rangle\langle\Psi,x\rangle\langle\Psi,y\rangle$ (360)
$\displaystyle=\sum_{x\in\\{0,1\\}^{4}_{+}}(-1)^{\vartheta(x)}\langle\Psi,x\rangle\langle\Psi,\overline{x}\rangle$
(361) $\displaystyle=\sum_{x\in
E}(-1)^{\vartheta(x)}\langle\Psi,x\rangle\langle\Psi,\overline{x}\rangle+(-1)^{\vartheta(\overline{x})}\langle\Psi,\overline{x}\rangle\langle\Psi,x\rangle$
(362) $\displaystyle=2\sum_{x\in
E}(-1)^{\vartheta(x)}\langle\Psi,x\rangle\langle\Psi,\overline{x}\rangle\ ,$
(363)

where the second step follows from $\theta u=\overline{u}\theta$ for all
$u\in\mathbb{C}$ due to the antiunitarity of $\theta$, the third step from Eq.
(204) and the fourth and final steps from
$\vartheta(\overline{x})=\vartheta(x)$ for $x\in\\{0,1\\}^{4}_{+}$. Thus,

$\displaystyle\langle\Psi,\theta\Psi\rangle=0\qquad\textrm{ if and only if
}\qquad\sum_{x\in\\{0,1\\}^{4}_{+}}(-1)^{\vartheta(x)}\langle\Psi,x\rangle\langle\Psi,\overline{x}\rangle=0\
.$ (364)

We proceed to prove that $\Lambda(\ket{\Psi}\otimes\ket{\Psi})=0$ is
equivalent to Eq. (364). We start by using Eq. (21) to write the operator
$\Lambda$ in terms of creation and annihilation operators:

$\displaystyle\Lambda$ $\displaystyle=\sum_{j=1}^{8}c_{j}\otimes c_{j}$ (365)
$\displaystyle=\sum_{j=1}^{4}\left(c_{2j-1}\otimes c_{2j-1}+c_{2j}\otimes
c_{2j}\right)$ (366)
$\displaystyle=\sum_{j=1}^{4}\left((a_{j}+a_{j}^{\dagger})\otimes(a_{j}+a_{j}^{\dagger})+i(a_{j}-a_{j}^{\dagger})\otimes
i(a_{j}-a_{j}^{\dagger})\right)$ (367)
$\displaystyle=2\sum_{j=1}^{4}\left(a_{j}\otimes
a_{j}^{\dagger}+a_{j}^{\dagger}\otimes a_{j}\right).$ (368)

Applying this expression to $\ket{\Psi}\otimes\ket{\Psi}$ and using Eq. (19)
gives

$\displaystyle\Lambda\left(\ket{\Psi}\otimes\ket{\Psi}\right)$
$\displaystyle=2\sum_{x,y\in\\{0,1\\}^{4}_{+}}\sum_{j=1}^{4}\left(a_{j}\otimes
a_{j}^{\dagger}+a_{j}^{\dagger}\otimes
a_{j}\right)\left(\ket{x}\otimes\ket{y}\right)\langle x,\Psi\rangle\langle
y,\Psi\rangle$ (369)
$\displaystyle=2\sum_{x,y\in\\{0,1\\}^{4}_{+}}\sum_{j=1}^{4}(-1)^{\eta_{j}(x+y)}(x_{j}\overline{y_{j}}+\overline{x_{j}}y_{j})\left(\ket{x\oplus
e_{j}}\otimes\ket{y\oplus e_{j}}\right)\langle x,\Psi\rangle\langle
y,\Psi\rangle$
$\displaystyle=2\sum_{x,y\in\\{0,1\\}^{4}_{+}}\sum_{j=1}^{4}(-1)^{\eta_{j}(x+y)}(y_{j}\oplus\overline{x_{j}})\left(\ket{x\oplus
e_{j}}\otimes\ket{y\oplus e_{j}}\right)\langle x,\Psi\rangle\langle
y,\Psi\rangle$
$\displaystyle=2\sum_{x,y\in\\{0,1\\}^{4}_{-}}\left(\sum_{j=1}^{4}(-1)^{\eta_{j}(x+y)}(y_{j}\oplus\overline{x_{j}})\langle
x\oplus e_{j},\Psi\rangle\langle y\oplus
e_{j},\Psi\rangle\right)\left(\ket{x}\otimes\ket{y}\right)$

where in the third line we used
$u\overline{v}+\overline{u}v=u\oplus\overline{v}$ for all $u,v\in\\{0,1\\}$
and in the last line we used $\left\\{x\oplus
e_{j}\,|\,x\in\\{0,1\\}^{4}_{+}\right\\}=\\{0,1\\}^{4}_{-}$ ,
$\eta_{j}(x\oplus e_{j})=\eta_{j}(x)$ and $u\oplus v\oplus 2e_{j}=u\oplus v$
valid for $x\in\\{0,1\\}_{+}^{4}$, $j\in[4]$ and $u,v\in\\{0,1\\}$. It follows
that $\Lambda\left(\ket{\Psi}\otimes\ket{\Psi}\right)=0$ if and only if

$\displaystyle\sum_{j=1}^{4}(-1)^{\eta_{j}(x+y)}(y_{j}\oplus\overline{x_{j}})\langle
x\oplus e_{j},\Psi\rangle\langle y\oplus e_{j},\Psi\rangle=0\quad\text{ for
all }\quad x,y\in\\{0,1\\}^{4}_{+}\ .$ (370)

Since $x$ and $y$ have the same parity, either $|x-y|=4$ (i.e.,
$y=\overline{x}$), $|x-y|=2$ or $|x-y|=0$ (i.e., $y=x$). The expression (370)
is non-zero only if $|x-y|=4$ (we argue below why). For $|x-y|=4$ Eq. (370)
becomes

$\displaystyle\sum_{j=1}^{4}(-1)^{\eta_{j}(x+\bar{x})}\langle x\oplus
e_{j},\Psi\rangle\langle\bar{x}\oplus e_{j},\Psi\rangle=0\quad\text{ if and
only if }\quad\sum_{z\in E}^{4}(-1)^{\vartheta(z)}\langle z,\Psi\rangle\langle
z,\Psi\rangle=0\ ,$ (371)

where we used $\\{x\oplus e_{j}\,|\,j\in[4]\\}=E$ with
$E\subset\\{0,1\\}_{+}^{4}$, $|E|=4$ a subset of even-weight strings such that
$E\cup\overline{E}=\\{0,1\\}^{4}_{+}$. We also used that
$(-1)^{\eta_{j}(x+\overline{x})}=(-1)^{j-1}$ can be replaced by
$(-1)^{\vartheta(x)+\vartheta(z)}$ upon changing the summation over $j\in[4]$
to a summation over $z\in E$. We recovered the right hand side of Eq. (364),
proving the claim.

It remains to argue that Eq. (370) is zero for $|x-y|\in\\{0,2\\}$. For
$|x-y|=0$, i.e., $x=y$, Eq. (370) is zero because
$x_{j}\otimes\overline{x_{j}}=0$ for $j\in[4]$. We exemplify that terms
$\propto\ket{x}\otimes\ket{y}$ with $|x-y|=2$ are zero by considering $x=1000$
and $y=0100$. Starting from Eq. (369) we obtain

$\displaystyle\left(\bra{1000}\otimes\bra{0100}\right)\Lambda\left(\ket{\Psi}\otimes\ket{\Psi}\right)$
$\displaystyle=((-1)^{\eta_{1}(1100)}+(-1)^{\eta_{2}(1100)})\langle
0000,\Psi\rangle\langle 1100,\Psi\rangle=0\ ,$ (372)

where $\eta_{1}(1100)=0$ and $\eta_{2}(1100)=1$. The remaining cases with
$|x-y|=2$ proceed similarly.

∎

##  Appendix B Commutativity of the map $\theta$ and quadratic Majorana
monomials

In the following, we prove Lemma 5.2.

###### Proof.

We start by showing that $\theta=c_{1}c_{3}c_{5}c_{7}K$, where $K$ denotes the
antiunitary given by complex conjugation in the number state basis. For this,
it suffices to show that $\theta$ is antiunitary, which directly follows from
unitarity of $c_{1}c_{3}c_{5}c_{7}$, and that it satisfies Eq. (204). We show
the later using Eqs. (21) and (19): We have

$\displaystyle\theta\ket{x}$ $\displaystyle=c_{1}c_{3}c_{5}c_{7}K\ket{x}$
(373)
$\displaystyle=(a_{1}+a_{1}^{\dagger})(a_{2}+a_{2}^{\dagger})(a_{3}+a_{3}^{\dagger})(a_{4}+a_{4}^{\dagger})\ket{x}$
(374)
$\displaystyle=(-1)^{\eta_{4}(x)+\eta_{3}(x)+\eta_{2}(x)+\eta_{1}(x)}(x_{1}+\overline{x_{1}})(x_{2}+\overline{x_{2}})(x_{3}+\overline{x_{3}})(x_{4}+\overline{x_{4}})\ket{\overline{x}}$
(375) $\displaystyle=(-1)^{\vartheta(x)}\ket{\overline{x}}\ ,$ (376)

where we used
$(-1)^{\eta_{4}(x)+\eta_{3}(x)+\eta_{2}(x)+\eta_{1}(x)}=(-1)^{3x_{1}+2x_{2}+x_{3}}=(-1)^{x_{1}+x_{3}}=(-1)^{\vartheta(x)}$
and $x_{j}+\overline{x_{j}}=1$ for $j\in[4]$.

The result $\theta c_{j}c_{k}=c_{j}c_{k}\theta~{}$ follows from simple algebra
considering $c_{2j}K=Kc_{2j}$ and $c_{2j-1}K=-Kc_{2j-1}$. We show these last
two equalities by explicitly computing their action on
$x\in\\{0,1\\}_{+}^{4}$:

$\displaystyle Kc_{2j}\ket{x}$
$\displaystyle=K(a_{j}+a_{j}^{\dagger})\ket{x}=(x_{j}+\overline{x_{j}})\ket{x\oplus
e_{j}}=\ket{x\oplus e_{j}}\ ,$ (377) $\displaystyle c_{2j}K\ket{x}$
$\displaystyle=(a_{j}+a_{j}^{\dagger})\ket{x}=(x_{j}+\overline{x_{j}})\ket{x\oplus
e_{j}}=Kc_{2j}\ket{x}$ (378)

and

$\displaystyle Kc_{2j+1}\ket{x}$
$\displaystyle=Ki(a_{j}-a_{j}^{\dagger})\ket{x}=-i(x_{j}-\overline{x_{j}})\ket{x\oplus
e_{j}}\ ,$ (379) $\displaystyle c_{2j+1}K\ket{x}$
$\displaystyle=i(a_{j}-a_{j}^{\dagger})\ket{x}=i(x_{j}-\overline{x_{j}})\ket{x\oplus
e_{j}}=-Kc_{2j+1}\ket{x}\ .$ (380)

∎

###### Proof.

We will prove that $\theta c_{j}=-c_{j}\theta$ for $j\in[8]$, which implies
the result.

We prove this for $j$ odd, the proof for $j$ even proceeds similarly. We use
Eq. (21) to write the Majorana operators as creation and annihilation
operators which act on basis states according to Eq. (19), and we apply
$\theta$ according to Eq. (204):

$\displaystyle\theta c_{2j-1}\ket{x}$
$\displaystyle=\theta(a_{j}+a_{j}^{\dagger})\ket{x}$ (381)
$\displaystyle=\theta(-1)^{\eta_{j}(x)}(x_{j}+\overline{x_{j}})\ket{x\oplus
e_{j}}$ (382) $\displaystyle=(-1)^{\eta_{j}(x)}(-1)^{\vartheta(x\oplus
e_{j})}(x_{j}+\overline{x_{j}})\ket{\overline{x\oplus e_{j}}}\ ,$ (383)
$\displaystyle c_{2j-1}\theta\ket{x}$
$\displaystyle=(-1)^{\vartheta(x)}(a_{j}+a_{j}^{\dagger})\ket{\overline{x}}$
(385)
$\displaystyle=(-1)^{\eta_{j}(\overline{x})}(-1)^{\vartheta(x)}(\overline{x_{j}}+x_{j})\ket{\overline{x}\oplus
e_{j}}\ .$ (386)

The equality $\theta c_{2j-1}=-c_{2j-1}\theta$ follows from
$\ket{\overline{x}\oplus e_{j}}=\ket{\overline{x\oplus e_{j}}}$ for $j\in[4]$,
from $(-1)^{\eta_{j}(\overline{x})}=(-1)^{j+1}(-1)^{\eta_{j}(x)}$ and from
$(-1)^{\vartheta(x)}=(-1)^{j}(-1)^{\vartheta(x\oplus e_{j})}$. ∎

## References

  * [1] Eric Dennis, Alexei Kitaev, Andrew Landahl, and John Preskill.  Topological quantum memory.  Journal of Mathematical Physics, 43(9):4452–4505, 2002.  https://doi.org/10.1063/1.1499754.
  * [2] Ruben S. Andrist, H. Bombin, Helmut G. Katzgraber, and M. A. Martin-Delgado.  Optimal error correction in topological subsystem codes.  Physical Review A, 85:050302, 2012.  https://doi.org/10.1103/PhysRevA.85.050302.
  * [3] Sergey Bravyi and Robert König.  Disorder-assisted error correction in Majorana chains.  Communications in Mathematical Physics, 316(3):641–692, 2012.  https://doi.org/10.1007/s00220-012-1606-9.
  * [4] Andrew S. Darmawan and David Poulin.  Tensor-network simulations of the surface code under realistic noise.  Physical Review Letters, 119:040502, 2017.  https://doi.org/10.1103/PhysRevLett.119.040502.
  * [5] Sergey Bravyi, Matthias Englbrecht, Robert König, and Nolan Peard.  Correcting coherent errors with surface codes.  npj Quantum Information, 4(1):1–6, 2018.  https://doi.org/10.1038/s41534-018-0106-y.
  * [6] David K. Tuckett, Andrew S. Darmawan, Christopher T. Chubb, Sergey Bravyi, Stephen D. Bartlett, and Steven T. Flammia.  Tailoring surface codes for highly biased noise.  Physical Review X, 9:041031, 2019.  https://doi.org/10.1103/PhysRevX.9.041031.
  * [7] Scott Aaronson and Daniel Gottesman.  Improved simulation of stabilizer circuits.  Physical Review A, 70(5):052328, 2004.  https://doi.org/10.1103/PhysRevA.70.052328.
  * [8] Leslie G. Valiant.  Quantum computers that can be simulated classically in polynomial time.  In Proceedings of the Thirty-Third Annual ACM Symposium on Theory of Computing, STOC ’01, page 114–123, New York, USA, 2001.  https://doi.org/10.1145/380752.380785.
  * [9] Barbara M. Terhal and David P. DiVincenzo.  Classical simulation of noninteracting-fermion quantum circuits.  Physical Review A, 65:032325, 2002.  https://doi.org/10.1103/PhysRevA.65.032325.
  * [10] Emanuel Knill.  Fermionic linear optics and matchgates.  Technical Report LAUR-01-4472, Los Alamos National Laboratory, 2001.
  * [11] Sergey Bravyi and Alexei Kitaev.  Universal quantum computation with ideal Clifford gates and noisy ancillas.  Physical Review A, 71:022316, 2005.  https://doi.org/10.1103/PhysRevA.71.022316.
  * [12] Sergey Bravyi.  Universal quantum computation with the $\nu=5/2$ fractional quantum Hall state.  Physical Review A, 73:042313, 2006.  https://doi.org/10.1103/PhysRevA.73.042313.
  * [13] Martin Hebenstreit, Richard Jozsa, Barbara Kraus, Sergii Strelchuk, and Mithuna Yoganathan.  All pure fermionic non-Gaussian states are magic states for matchgate computations.  Physical Review Letters, 123:080503, 2019.  https://doi.org/10.1103/PhysRevLett.123.080503.
  * [14] Mark Howard and Earl Campbell.  Application of a resource theory for magic states to fault-tolerant quantum computing.  Physical Review Letters, 118:090501, 2017.  https://doi.org/10.1103/PhysRevLett.118.090501.
  * [15] Markus Heinrich and David Gross.  Robustness of magic and symmetries of the stabiliser polytope.  Quantum, 3:132, 2019.  https://doi.org/10.22331/q-2019-04-08-132.
  * [16] Sergey Bravyi, Graeme Smith, and John A. Smolin.  Trading classical and quantum computational resources.  Physical Review X, 6(2):021043, 2016.  https://doi.org/10.1103/PhysRevX.6.021043.
  * [17] Sergey Bravyi and David Gosset.  Improved classical simulation of quantum circuits dominated by Clifford gates.  Physical Review Letters, 116(25), 2016.  https://doi.org/10.1103/PhysRevLett.116.250501.
  * [18] Sergey Bravyi, Dan Browne, Padraic Calpin, Earl Campbell, David Gosset, and Mark Howard.  Simulation of quantum circuits by low-rank stabilizer decompositions.  Quantum, 3:181, 2019.  https://doi.org/10.22331/q-2019-09-02-181.
  * [19] Arne Heimendahl, Felipe Montealegre-Mora, Frank Vallentin, and David Gross.  Stabilizer extent is not multiplicative.  Quantum, 5:400, 2021.  https://doi.org/10.22331/q-2021-02-24-400.
  * [20] Michael Beverland, Earl Campbell, Mark Howard, and Vadym Kliuchnikov.  Lower bounds on the non-Clifford resources for quantum computations.  Quantum Science and Technology, 5(3):035009, 2020.  https://doi.org/10.1088/2058-9565/ab8963.
  * [21] James R. Seddon, Bartosz Regula, Hakop Pashayan, Yingkai Ouyang, and Earl T. Campbell.  Quantifying quantum speedups: Improved classical simulation from tighter magic monotones.  PRX Quantum, 2:010345, 2021.  https://doi.org/10.1103/PRXQuantum.2.010345.
  * [22] Kaifeng Bu, Weichen Gu, and Arthur Jaffe.  Stabilizer testing and magic entropy.  arXiv:2306.09292, 2023.
  * [23] Joel A. Tropp.  Greed is good: algorithmic results for sparse approximation.  IEEE Transactions on Information Theory, 50(10):2231–2242, 2004\.  https://doi.org/10.1109/TIT.2004.834793.
  * [24] Scott Shaobing Chen, David L. Donoho, and Michael A. Saunders.  Atomic decomposition by basis pursuit.  SIAM Review, 43(1):129–159, 2001.  https://doi.org/10.1137/S003614450037906X.
  * [25] Jean-Jacques Fuchs.  On sparse representations in arbitrary redundant bases.  IEEE Transactions on Information Theory, 50(6):1341–1344, 2004\.  https://doi.org/10.1109/TIT.2004.828141.
  * [26] Joel A. Tropp.  Recovery of short, complex linear combinations via /spl lscr//sub 1/ minimization.  IEEE Transactions on Information Theory, 51(4):1568–1570, 2005\.  https://doi.org/10.1109/TIT.2005.844057.
  * [27] Hakop Pashayan, Oliver Reardon-Smith, Kamil Korzekwa, and Stephen D. Bartlett.  Fast estimation of outcome probabilities for quantum circuits.  PRX Quantum, 3:020361, 2022.  https://doi.org/10.1103/PRXQuantum.3.020361.
  * [28] Farid Alizadeh and Donald Goldfarb.  Second-order cone programming.  Mathematical Programming, 95(1):3–51, 2003.  https://doi.org/10.1007/s10107-002-0339-5.
  * [29] Stephen P. Boyd and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.  https://doi.org/10.1017/CBO9780511804441.
  * [30] Venkat Chandrasekaran, Benjamin Recht, Pablo A. Parrilo, and Alan S. Willsky.  The Convex Geometry of Linear Inverse Problems.  Foundations of Computational Mathematics, 12(6):805–849, 2012.  https://doi.org/10.1007/s10208-012-9135-7.
  * [31] Christophe Piveteau and David Sutter.  Circuit knitting with classical communication.  IEEE Transactions on Information Theory, 70(4):2734–2745, 2024\.  https://doi.org/10.1109/TIT.2023.3310797.
  * [32] Fernando de Melo, Piotr Ćwikliński, and Barbara M. Terhal.  The power of noisy fermionic quantum computation.  New Journal of Physics, 15(1):013015, 2013.  https://doi.org/10.1088/1367-2630/15/1/013015.
  * [33] Richard Jozsa and Akimasa Miyake.  Matchgates and classical simulation of quantum circuits.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 464(2100):3089–3106, 2008.  https://doi.org/10.1098/rspa.2008.0189.
  * [34] Stephen D. Bartlett and Barry C. Sanders.  Efficient classical simulation of optical quantum information circuits.  Physical Review Letters, 89:207903, 2002.  https://doi.org/10.1103/PhysRevLett.89.207903.
  * [35] Stephen D. Bartlett, Barry C. Sanders, Samuel L. Braunstein, and Kae Nemoto.  Efficient classical simulation of continuous variable quantum information processes.  Physical Review Letters, 88:097904, 2002.  https://doi.org/10.1103/PhysRevLett.88.097904.
  * [36] Rolando Somma, Howard Barnum, Gerardo Ortiz, and Emanuel Knill.  Efficient solvability of Hamiltonians and limits on the power of some quantum computational models.  Physical Review Letters, 97:190501, 2006.  https://doi.org/10.1103/PhysRevLett.97.190501.
  * [37] Daniel Gottesman.  Stabilizer Codes and Quantum Error Correction.  PhD thesis, California Institute of Technology, 1997.  https://doi.org/10.7907/rzr7-dt72.
  * [38] Victor Veitch, Christopher Ferrie, David Gross, and Joseph Emerson.  Negative quasi-probability as a resource for quantum computation.  New Journal of Physics, 14(11):113011, 2012.  https://doi.org/10.1088/1367-2630/14/11/113011.
  * [39] Andrea Mari and Jens Eisert.  Positive wigner functions render classical simulation of quantum computation efficient.  Physical Review Letters, 109:230503, 2012.  https://doi.org/10.1103/PhysRevLett.109.230503.
  * [40] Hakop Pashayan, Joel J. Wallman, and Stephen D. Bartlett.  Estimating outcome probabilities of quantum circuits using quasiprobabilities.  Physical Review Letters, 115:070501, 2015.  https://doi.org/10.1103/PhysRevLett.115.070501.
  * [41] Robert Raussendorf, Juani Bermejo-Vega, Emily Tyhurst, Cihan Okay, and Michael Zurel.  Phase-space-simulation method for quantum computation with magic states on qubits.  Physical Review A, 101:012350, 2020.  https://doi.org/10.1103/PhysRevA.101.012350.
  * [42] Robert Raussendorf, Dan E. Browne, Nicolas Delfosse, Cihan Okay, and Juan Bermejo-Vega.  Contextuality and wigner-function negativity in qubit quantum computation.  Physical Review A, 95:052334, 2017.  https://doi.org/10.1103/PhysRevA.95.052334.
  * [43] Juan Bermejo-Vega, Nicolas Delfosse, Dan E. Browne, Cihan Okay, and Robert Raussendorf.  Contextuality as a resource for models of quantum computation with qubits.  Physical Review Letters, 119:120505, 2017.  https://doi.org/10.1103/PhysRevLett.119.120505.
  * [44] Markus Frembs, Sam Roberts, and Stephen D. Bartlett.  Contextuality as a resource for measurement-based quantum computation beyond qubits.  New Journal of Physics, 20(10):103011, 2018.  https://doi.org/10.1088/1367-2630/aae3ad.
  * [45] Anna Vershynina.  Complete criterion for convex-gaussian-state detection.  Physical Review A, 90:062329, 2014.  https://doi.org/10.1103/PhysRevA.90.062329.
  * [46] Shigeo Hakkaku, Yuichiro Tashima, Kosuke Mitarai, Wataru Mizukami, and Keisuke Fujii.  Quantifying fermionic nonlinearity of quantum circuits.  Physical Review Research, 4:043100, 2022.  https://doi.org/10.1103/PhysRevResearch.4.043100.
  * [47] Avinash Mocherla, Lingling Lao, and Dan E. Browne.  Extending matchgate simulation methods to universal quantum circuits.  arxiv:2302.02654, 2023.
  * [48] Oliver Reardon-Smith, Michał Oszmaniec, and Kamil Korzekwa.  Improved classical simulation of quantum circuits dominated by fermionic linear optical gates.  arXiv:2307.12702, 2023.
  * [49] Joshua Cudby and Sergii Strelchuk.  Gaussian decomposition of magic states for matchgate computations.  arXiv:2307.12654, 2023.
  * [50] Alonso Botero and Benni Reznik.  BCS-like modewise entanglement of fermion Gaussian states.  Physics Letters A, 331(1):39–44, 2004.  https://doi.org/10.1016/j.physleta.2004.08.037.
  * [51] Grigori G. Amosov and Sergey N. Filippov.  Spectral properties of reduced fermionic density operators and parity superselection rule.  Quantum Information Processing, 16(1):2, 2016.  https://doi.org/10.1007/s11128-016-1467-9.
  * [52] Wallace Givens.  Computation of plain unitary rotations transforming a general matrix to triangular form.  Journal of the Society for Industrial and Applied Mathematics, 6(1):26–50, 1958.  https://doi.org/10.1137/0106004.
  * [53] Gene H. Golub and Charles F. Van Loan.  Matrix Computations (3rd Ed.).  Johns Hopkins University Press, USA, 1996.  https://doi.org/10.56021/9781421407944.
  * [54] Sergey Bravyi.  Lagrangian representation for fermionic linear optics.  Quantum Information & Computation, 5(3):216–238, 2005.  https://doi.org/10.26421/qic5.3-3.
  * [55] Woody Lichtenstein.  A system of quadrics describing the orbit of the highest weight vector.  Proceedings of the American Mathematical Society, 84(4):605–608, 1982.  https://doi.org/10.1090/s0002-9939-1982-0643758-8.
  * [56] Marek Kuś and Ingemar Bengtsson.  “Classical” quantum states.  Physical Review A, 80(2):022319, 2009.  https://doi.org/10.1103/PhysRevA.80.022319.
  * [57] Michał Oszmaniec and Marek Kuś.  On detection of quasiclassical states.  Journal of Physics A: Mathematical and Theoretical, 45(24):244034, 2012.  https://doi.org/10.1088/1751-8113/45/24/244034.
  * [58] Michał Oszmaniec, Jan Gutt, and Marek Kuś.  Classical simulation of fermionic linear optics augmented with noisy ancillas.  Physical Review A, 90(2):020302, 2014.  https://doi.org/10.1103/PhysRevA.90.020302.
  * [59] Per-Olov Löwdin.  Quantum theory of many-particle systems. i. physical interpretations by means of density matrices, natural spin-orbitals, and convergence problems in the method of configurational interaction.  Physical Review, 97:1474–1489, 1955.  https://doi.org/10.1103/PhysRev.97.1474.
  * [60] Sergey Bravyi and David Gosset.  Complexity of quantum impurity problems.  Communications in Mathematical Physics, 356(2), 2017.  https://doi.org/10.1007/s00220-017-2976-9.
  * [61] Huzihiro Araki and Hajime Moriya.  Joint extension of states of subsystems for a CAR system.  Communications in Mathematical Physics, 237(1):105–122, 2003.  https://doi.org/10.1007/s00220-003-0832-6.
  * [62] Thomas P Hayes.  A large-deviation inequality for vector-valued martingales, 2005.  Unpublished manuscript.  https://www.cs.unm.edu/$\sim$hayes/papers/VectorAzuma/.
  * [63] Gongguo Tang, Badri Narayan Bhaskar, and Benjamin Recht.  Sparse recovery over continuous dictionaries-just discretize.  In 2013 Asilomar Conference on Signals, Systems and Computers, pages 1043–1047, 2013.  https://doi.org/10.1109/ACSSC.2013.6810450.