# Classical simulation of non-Gaussian fermionic circuits

Beatriz Dias  Department of Mathematics, School of Computation, Information
and Technology, Technical University of Munich, Garching, Germany  Munich
Center for Quantum Science and Technology, Munich, Germany  Robert König
Department of Mathematics, School of Computation, Information and Technology,
Technical University of Munich, Garching, Germany  Munich Center for Quantum
Science and Technology, Munich, Germany

###### Abstract

We propose efficient algorithms for classically simulating fermionic linear
optics operations applied to non-Gaussian initial states. By gadget
constructions, this provides algorithms for fermionic linear optics with non-
Gaussian operations. We argue that this problem is analogous to that of
simulating Clifford circuits with non-stabilizer initial states: Algorithms
for the latter problem immediately translate to the fermionic setting. Our
construction is based on an extension of the covariance matrix formalism which
permits to efficiently track relative phases in superpositions of Gaussian
states. It yields simulation algorithms with polynomial complexity in the
number of fermions, the desired accuracy, and certain quantities capturing the
degree of non-Gaussianity of the initial state. We study one such quantity,
the fermionic Gaussian extent, and show that it is multiplicative on tensor
products when the so-called fermionic Gaussian fidelity is. We establish this
property for the tensor product of two arbitrary pure states of four fermions
with positive parity.

###### Contents

  1. 1 Introduction
    1. 1.1 Efficiently simulable quantum computations
      1. 1.1.1 Clifford circuits / Stabilizer computations
      2. 1.1.2 Fermionic linear optics / Fermionic Gaussian computations
    2. 1.2 Classical simulation algorithms and measures of magic
    3. 1.3 The fermionic Gaussian extent
    4. 1.4 On the (sub)multiplicativity of the extent
    5. 1.5 Our contribution
    6. 1.6 Prior and related work
    7. 1.7 Outline
  2. 2 Background
    1. 2.1 Dirac and Majorana operators
    2. 2.2 Gaussian unitaries
    3. 2.3 Fermionic Gaussian (pure) states
    4. 2.4 Gaussianity condition
    5. 2.5 Covariance matrices, Gaussian states and Wick’s theorem
    6. 2.6 Inner product formulas for Gaussian states
    7. 2.7 Gaussian evolution and occupation number measurement
    8. 2.8 The tensor product of two fermionic states
  3. 3 Tracking relative phases in fermionic linear optics
    1. 3.1 Subroutines
    2. 3.2 Computing overlaps and descriptions of evolved/measured states
    3. 3.3 Initial states for computation
  4. 4 Classical simulation of fermionic Gaussian circuits with non-Gaussian initial states
    1. 4.1 Extending simulation algorithms to superpositions
    2. 4.2 Sparsification: Relating $\mathcal{D}$-extent to approximate $\mathcal{D}$-rank
    3. 4.3 Fast norm estimation and approximate simulation
    4. 4.4 Fermionic linear optics with non-Gaussian initial states
    5. 4.5 Efficient additive-error strong simulation
  5. 5 Multiplicativity of the Gaussian fidelity for $4$ fermions
    1. 5.1 Four-fermion Gaussian and non-Gaussian states
    2. 5.2 The Gaussian fidelity for 4-fermion states
    3. 5.3 Multiplicativity of the Gaussian fidelity for 4-fermion states
  6. 6 Multiplicativity of $\mathcal{D}$-fidelity implies that of $\mathcal{D}$-extent
    1. 6.1 Multiplicativity for finite dictionaries
    2. 6.2 Multiplicativity for infinite dictionaries
  7. 7 Multiplicativity of the Gaussian extent for four fermions
  8. A Alternative Gaussianity condition for $4$-fermion states
  9. B Commutativity of the map $\theta$ and quadratic Majorana monomials

##  1 Introduction

While universal polynomial-time quantum computation is believed to exceed the
capabilities of efficient classical algorithms, restricted classes of quantum
computations are amenable to efficient classical simulation. Identifying such
models and corresponding simulation algorithms is a central goal in the study
of quantum computing. On the one hand, a good characterization of the boundary
between the computational power of classical and quantum computational models
provides insight into potential quantum advantages. On the other hand,
efficient classical simulation methods can be used to assess the merits and
scalability of quantum information-processing proposals. For example, the
resilience of certain quantum codes against restricted noise models has
successfully been studied by means of classical simulation methods giving
threshold estimates for large-scale systems, see e.g., [1, 2, 3, 4, 5, 6] for
an incomplete list of relevant references.

###  1.1 Efficiently simulable quantum computations

Most known examples of efficiently simulable quantum computations can be
summarized by the following ingredients:

  1. (i)

A set $\mathcal{D}$ of states with the property that each element
$\Psi\in\mathcal{D}$ has a succinct classical description $d_{\Psi}$. In the
following, we will refer to $\mathcal{D}$ as a dictionary.

  2. (ii)

A set $\mathcal{E}$ of operations (unitary or non-unitary evolutions), again
with the property that each element $E\in\mathcal{E}$ has a succinct classical
description $d_{E}$. Following resource-theoretic conventions, we call
$\mathcal{E}$ the set of free operations.

  3. (iii)

A set $\mathcal{M}$ of measurements (quantum instruments) with an efficient
classical descriptions $d_{M}$ for each $M\in\mathcal{M}$, and the property
that every post-measurement state (associated with different measurement
outcomes) obtained by applying $M\in\mathcal{M}$ to a state
$\Psi\in\mathcal{D}$ belongs to $\mathcal{D}$.

A triple $(\mathcal{D},\mathcal{E},\mathcal{M})$ gives rise to a (generally
restricted) quantum computational model by composing these ingredients. A
typical (non-adaptive) computation proceeds by preparing an initial state
$\Psi\in\mathcal{D}$, applying a sequence
$\\{E_{t}\\}_{t=1}^{T}\subset\mathcal{E}$ of operations, and performing
measurements $\\{M_{k}\\}_{k=1}^{L}\subset\mathcal{M}$ in succession. Assuming
for simplicity that $\mathcal{E}$ consists of a set of unitaries, and that for
each $k\in[L]=\\{1,\dots,L\\}$, the measurement $M_{k}$ realizes a POVM
$M_{k}=\\{M^{(k)}_{m}\\}_{m\in\mathcal{M}_{k}}$ with outcomes in a set
$\mathcal{M}_{k}$, such a computation produces a sample
$m=(m_{1},\ldots,m_{L})\in\mathcal{M}_{1}\times\cdots\times\mathcal{M}_{L}$
from the distribution

$\displaystyle p(m_{1},\ldots,m_{L})$
$\displaystyle=\langle\Psi_{T},M^{(L)}_{m_{L}}\cdots
M^{(1)}_{m_{1}}\Psi_{T}\rangle\qquad\textrm{ where
}\qquad\Psi_{T}=E_{T}\circ\cdots\circ E_{1}(\Psi)\ .$ (1)

More generally, one may consider circuits where operations are chosen
adaptively depending on intermediate measurement results, assuming that the
dependence is given by an efficiently computable function.

The task of classically simulating the computational model associated with
$(\mathcal{D},\mathcal{E},\mathcal{M})$ comes in two flavors. The input in
both cases is the collection
$(d_{\Psi},\\{d_{E_{t}}\\}_{t=1}^{T},\\{d_{M_{k}}\\}_{k=1}^{L})$ of
descriptions of the initial state, the set of operations applied, and the
measurements. The problem of weak simulation then consists in producing a
sample $m\in\mathcal{M}_{1}\times\cdots\times\mathcal{M}_{L}$ drawn according
to the (ideal) output distribution $p(m)$ of the circuit given by Eq. (1). In
contrast, the problem of strong simulation consists of computing the output
probability $p(m)$ for a given (potential) measurement outcome $m$.

Relaxing the requirements of weak and strong simulation, one may allow for an
approximation error. For weak simulation, this is typically formalized by
demanding that the (probabilistic) classical algorithm outputs a sample $m$
drawn from a distribution $\tilde{p}$ which is $\delta$-close in
$L^{1}$-distance (for a chosen error parameter $\delta>0$) to the ideal output
distribution $p$. Similarly, for strong simulation, the output $\tilde{p}$ is
required to be close to the value $p(m)$ with a controlled (additive or
multiplicative) error.

In cases where a computational model specified by
$(\mathcal{D},\mathcal{E},\mathcal{M})$ is amenable to efficient classical
simulation, associated classical simulation algorithms are typically
constructed by considering evolution and measurement separately. The basic
problem then consists in constructing efficient classical algorithms with the
following functionalities:

  1. (a)

An algorithm $\mathsf{evolve}$ which, given classical descriptions $d_{\Psi}$
of a state $\Psi\in\mathcal{D}$ and $d_{E}$ of an evolution operation
$E\in\mathcal{E}$, computes a classical description $d_{E(\Psi)}$ of the
evolved state $E(\Psi)$.

  2. (b)

Given a classical description $d_{\Psi}$ of a state $\Psi\in\mathcal{D}$, a
classical description $d_{M}$ of a measurement $M\in\mathcal{M}$ (with
associated set of measurement outcomes $\mathcal{M}_{M}$) and a measurement
outcome $m\in\mathcal{M}_{M}$,

    1. (b1)

an algorithm $\mathsf{measureprob}$ which outputs the probability $p(m)$
(determined by Born’s rule) of obtaining measurement outcome $m$, and

    2. (b2)

an algorithm $\mathsf{postmeasure}$ which outputs a classical description of
the post-measurement state

when applying the measurement $M$ to $\Psi$.

It is clear that a triple
$(\mathsf{evolve},\mathsf{measureprob},\mathsf{postmeasure})$ of such
algorithms immediately gives rise to an efficient algorithm for strong
simulation of the model $(\mathcal{D},\mathcal{E},\mathcal{M})$, with a
runtime

$\displaystyle
T\cdot\mathsf{time}(\mathsf{evolve})+L\cdot(\mathsf{time}(\mathsf{measureprob})+\mathsf{time}(\mathsf{postmeasure}))$
(2)

which is linear in the number $T$ of operations applied, and linear in the
number $L$ of measurements. Assuming that for any measurement
$M\in\mathcal{M}$, the set of measurement outcomes $\mathcal{M}_{M}$
associated with $M$ is of constant cardinality, the triple
$(\mathsf{evolve},\mathsf{measureprob},\allowbreak\mathsf{postmeasure})$ also
gives rise to a randomized algorithm for weak simulation: Such an algorithm is
obtained by using $\mathsf{measureprob}$ to compute the entire distribution
$\\{p(m)\\}_{m\in\mathcal{M}_{M}}$ of measurement outcomes (when applying a
measurement $M$), and then drawing $m\in\mathcal{M}_{M}$ randomly according to
this distribution. The runtime of this probabilistic algorithm is

$\displaystyle
T\cdot\mathsf{time}(\mathsf{evolve})+L\cdot(\mathsf{time}(\mathsf{measureprob})\cdot
w+\mathsf{time}(\mathsf{postmeasure}))$ (3)

where $w=\max_{M\in\mathcal{M}}|\mathcal{M}_{M}|$ bounds the maximal
cardinality of the set of measurement outcomes.

####  1.1.1 Clifford circuits / Stabilizer computations

Perhaps the most well-known example of a computational model
$(\mathcal{D},\mathcal{E},\mathcal{M})$ where efficient algorithms
$(\mathsf{evolve},\mathsf{measureprob},\mathsf{postmeasure})$ can be provided
is the Gottesman-Knill-theorem for stabilizer computations on $n$ qubits. Here
$\mathcal{D}$ is the set $\mathsf{STAB}_{n}$ of $n$-qubit stabilizer states
(whose elements can be specified by their stabilizer generators, i.e.,
corresponding stabilizer tableaux), $\mathcal{E}$ is the set of Clifford
unitaries (described by symplectic matrices), and $\mathcal{M}$ are
measurements of single-qubit Pauli $Z$ operators (described by an index
$j\in[n]$). In this case, there are efficient algorithms with runtimes given
in Table 1.

algorithm | $\mathsf{time}$
---|---
$\mathsf{evolve}$ | $O(n)$
$\mathsf{measureprob}$ | $O(n)$
$\mathsf{postmeasure}$ | $O(n^{2})$
Table 1: Runtimes of building blocks $\mathsf{evolve}$,
$\mathsf{measureprob}$, $\mathsf{postmeasure}$ for classical simulation of
$n$-qubit stabilizer circuits as given in [7]. Evolution corresponds to
application of an $n$-qubit Clifford unitary, and each measurement is that of
a Pauli observable $Z_{j}$ with $j\in[n]$.

####  1.1.2 Fermionic linear optics / Fermionic Gaussian computations

A different class of efficiently simulable computations – the one we are
interested in here – is that of fermionic linear optics on $n$ fermions. We
focus on pure-state computations: Here the dictionary $\mathcal{D}$ consists
of the set $\mathcal{G}_{n}$ of pure fermionic Gaussian states. An element
$\Psi\in\mathcal{G}_{n}$ in the dictionary can be described by its covariance
matrix $\Gamma_{\Psi}$, an antisymmetric $2n\times 2n$ matrix with real
entries. The set $\mathcal{E}=\mathcal{E}_{\textrm{Gauss}}$ can be taken as
the set of Gaussian unitary operations. Each such unitary $U=U_{R}$ is fully
determined by an element $R\in O(2n)$ of the orthogonal group on
$\mathbb{R}^{2n}$, where $R\mapsto U_{R}$ defines a (projective) unitary
representation of $O(2n)$ on the space $\mathcal{H}^{n}$ of $n$ fermions. The
set $\mathcal{M}=\mathcal{M}_{\textrm{number}}$ consists of all occupation
number measurements. As in the case of stabilizer states, there are
polynomial-time algorithms
$(\mathsf{evolve},\mathsf{measureprob},\mathsf{postmeasure})$ for classical
simulation with runtimes summarized in Table 2. In particular, the covariance
matrix $\Gamma_{U_{R}\Psi}$ of a Gaussian state $\Psi$ evolved under a
Gaussian unitary $U_{R}$ can be computed in time $O(n^{3})$ from
$\Gamma_{\Psi}$ and $R$. The outcome probability of observing $0$
(respectively $1$) when performing an occupation number measurement can be
computed in time $O(1)$, and the covariance matrix of the post-measurement
state can be computed in time $O(n^{2})$ [8, 9, 10] (see also [3]).

algorithm | $\mathsf{time}$
---|---
$\mathsf{evolve}$ | $O(n^{3})$
$\mathsf{measureprob}$ | $O(1)$
$\mathsf{postmeasure}$ | $O(n^{2})$
Table 2: Runtimes of building blocks $\mathsf{evolve}$,
$\mathsf{measureprob}$, $\mathsf{postmeasure}$ for classical simulation of
$n$-fermion linear optics circuits as proposed in [8, 9, 10], see also [3].
Evolution amounts to application of a fermionic Gaussian unitary. Measurement
corresponds to measuring an observable $a_{j}^{\dagger}a_{j}$ (occupation
number) for $j\in[n]$.

###  1.2 Classical simulation algorithms and measures of magic

A natural way of extending the power of a quantum computational model
specified by $(\mathcal{D},\mathcal{E},\mathcal{M})$ consists in providing
resources/capabilities that do not belong to the specified sets. “Magic
states” are a prime example: Here a state $\Psi\not\in\mathcal{D}$ not
belonging to the dictionary is provided as an (initial) state in the quantum
computation, thereby providing additional capabilities to the computational
model. For example, non-Clifford unitaries can be realized by certain
stabilizer-computations (sometimes referred to as “gadgets”) applied to so-
called magic states [11]. Similarly, non-Gaussian initial states can be
combined with fermionic linear optics operations to realize non-Gaussian
operations [12, 13]. While such a magic state can even promote the
computational model to universal quantum computation, this is generally not
the case for all states $\Psi$. It is thus a natural question to quantify the
degree of “magicness” provided by a state $\Psi\not\in\mathcal{D}$. For the
set $\mathsf{STAB}_{n}$ of $n$-qubit stabilizer states, corresponding magic
monotones considered in the literature include the robustness of magic [14,
15], the exact and approximate stabilizer rank [16, 17, 18], the stabilizer
extent [18, 19], the stabilizer nullity [20], the generalized robustness [21]
and the magic entropy [22].

The maximum overlap of a given state $\Psi$ with an element of the dictionary
$\mathcal{D}$, i.e., the quantity

$\displaystyle F_{\mathcal{D}}(\Psi)$
$\displaystyle=\sup_{\varphi\in\mathcal{D}}|\langle\varphi,\Psi\rangle|^{2}\
,$ (4)

is arguably one of the most direct ways of quantifying how far $\Psi$ is from
a “free” state, i.e., a state belonging to $\mathcal{D}$. Motivated by the
analogously defined notion of stabilizer fidelity in Ref. [18], we call
$F_{\mathcal{D}}(\Psi)$ the $\mathcal{D}$-fidelity of $\Psi$ in the following.
This quantity plays an important role in our arguments when considering
multiplicativity properties. However, the $\mathcal{D}$-fidelity
$F_{\mathcal{D}}(\Psi)$ is not a good quantifier of hardness of classical
simulation because simply replacing $\Psi$ by an element of $\mathcal{D}$
typically leads to a significant approximation error.

From the point of view of classical simulation, a relevant magicness measure
for a state $\Psi\not\in\mathcal{D}$ relates to the (added) complexity when
trying to simulate a quantum computation with initial state $\Psi$, built from
a triple $(\mathcal{D},\mathcal{E},\mathcal{M})$ allowing for efficient
classical simulation. One such measure, introduced in Ref. [16] for the case
of stabilizer computations, is the $\mathcal{D}$-rank
$\chi_{\mathcal{D}}(\Psi)$ of $\Psi$. (For $\mathcal{D}=\mathsf{STAB}_{n}$,
this is called the stabilizer rank of $\Psi$.) It is defined as the minimum
number of terms when decomposing $\Psi$ as a linear combination of elements of
$\mathcal{D}$, i.e.,

$\displaystyle\chi_{\mathcal{D}}(\Psi)=\min\left\\{\chi\in\mathbb{N}\ |\
\exists\\{\varphi_{j}\\}_{j=1}^{\chi}\subset\mathcal{D},\\{\gamma_{j}\\}_{j=1}^{\chi}\subset\mathbb{C}\textrm{
such that }\Psi=\sum_{j=1}^{\chi}\gamma_{j}\varphi_{j}\ \right\\}\ .$ (5)

In the context of signal processing, the corresponding optimization problem is
referred to as a sparse approximation problem. The $\mathcal{D}$-rank
$\chi_{\mathcal{D}}(\Psi)$ appears naturally when constructing and analyzing
simulation algorithms, but it suffers from a number of shortcomings: On the
one hand, the set of states $\Psi\in\mathcal{H}$ whose $\mathcal{D}$-rank is
less than the dimension of the Hilbert space $\mathcal{H}$ is a set of zero
Lebesgue measure [23, Proposition 4.1]. On the other hand, the quantity
$\chi_{\mathcal{D}}(\Psi)$ relates to the classical simulation cost of exactly
simulating dynamics involving the state $\Psi$. In practice, some
approximation error is typically acceptable, and corresponding simulations can
be achieved with lower cost. In other words, the quantity
$\chi_{\mathcal{D}}(\Psi)$ does not accurately reflect the cost of approximate
simulation.

A more operationally relevant quantity is the $\delta$-approximate
$\mathcal{D}$-rank $\chi_{\mathcal{D}}^{\delta}(\Psi)$ of $\Psi$ introduced in
Ref. [17], again for stabilizer computations. For a fixed approximation error
$\delta>0$, this is given by the minimum $\mathcal{D}$-rank of any state
$\Psi^{\prime}$ that is $\delta$-close to $\Psi$, i.e.,

$\displaystyle\chi_{\mathcal{D}}^{\delta}(\Psi)=\min\left\\{\chi_{\mathcal{D}}(\Psi^{\prime})\
|\ \Psi^{\prime}\in\mathcal{H}\text{ such that
}\|\Psi-\Psi^{\prime}\|\leq\delta\right\\}\ .$ (6)

An exact classical simulation algorithm whose complexity scales with the exact
$\mathcal{D}$-rank $\chi_{\mathcal{D}}(\Psi)$ provides an approximate
simulation at a cost with an identical scaling in the approximate (instead of
exact) $\mathcal{D}$-rank $\chi_{\mathcal{D}}^{\delta}(\Psi)$ of $\Psi$. Here
approximate weak simulation means that instead of sampling from the ideal
output distribution $P$ of a circuit, the simulation samples from a
distribution $P^{\prime}$ whose $L^{1}$-distance from $P$ is bounded by
$O(\delta)$. Similarly, in approximate (strong) simulation, output
probabilities are approximately computed with a controlled approximation
error.

A different quantity of interest is obtained by replacing the ill-behaved rank
function (i.e., size of the support) in the definition of the
$\mathcal{D}$-rank $\chi_{\mathcal{D}}(\Psi)$ by the $L^{1}$-norm of the
coefficients when representing $\Psi$ as a linear combination. In the context
of stabilizer states the corresponding quantity was introduced by Bravyi et
al. [18] under the term stabilizer extent: For a state
$\Psi\in(\mathbb{C}^{2})^{\otimes n}$ it is defined as

$\displaystyle\xi_{\mathsf{STAB}_{n}}(\Psi)$
$\displaystyle=\inf\left\\{\|\gamma\|_{1}^{2}\ |\
\gamma:\mathsf{STAB}_{n}\rightarrow\mathbb{C}\textrm{ such that
}\Psi=\sum_{\varphi\in\mathsf{STAB}_{n}}\gamma(\varphi)\varphi\right\\}\ ,$
(7)

where $\|\gamma\|_{1}=\sum_{\varphi\in\mathsf{STAB}_{n}}|\gamma(\varphi)|$
denotes the $1$-norm of $\gamma$. The corresponding convex optimization
problem is known as the basis pursuit problem [24] (when $\mathsf{STAB}_{n}$
is replaced by e.g., a finite dictionary $\mathcal{D}$). Sufficient conditions
for when the basis pursuit problem yields a solution of the sparse
approximation problem where investigated in a series of works culminating in
Fuchs’ condition [25] (see also [26]). More importantly for (approximate)
simulation, feasible solutions of the basis pursuit problem provide upper
bounds on the sparse approximation problem. For the stabilizer rank, a
sparsification result (see [18, Theorem 1]) gives an upper bound on the
$\delta$-approximate stabilizer rank $\chi_{\mathsf{STAB}_{n}}(\Psi)$ in terms
of the stabilizer extent $\xi_{\mathsf{STAB}_{n}}(\Psi)$, for any $\delta>0$
(see Section 4.2).

Building on earlier results [17], it was shown in Ref. [18] that a stabilizer
circuit on $n$ qubits with $L$ Clifford gates initialized in a state $\Psi$
can be weakly simulated with error $\delta$ in a time scaling as
$O({{\xi_{\mathsf{STAB}_{n}}(\Psi)}/{\delta^{2}}}\cdot\mathsf{poly}(n,L))$.
The error $\delta$ expresses the $L^{1}$-norm distance of the distribution of
simulated measurement outcomes from the output distribution of the actual
quantum computation. Here we are not accounting for the time required to
perform classical computations when adaptive quantum circuits are considered.
In addition, Ref. [27] provided a classical algorithm for strong simulation of
a circuit $U$ with $L$ Clifford gates and $t$ $T$-gates initialized in a
stabilizer state $\Psi$ with an additive error $\delta$. Their algorithm
outputs an estimate of the probability $|\langle x,U\Psi\rangle|^{2}$ of
obtaining measurement outcome $x\in\\{0,1\\}^{n}$ up to an additive error
$\delta$, with success probability greater than $1-p_{f}$. It has runtime
$O(\xi_{\mathsf{STAB}_{n}}(\ket{T}^{\otimes
t})\log(1/p_{f})\cdot\mathsf{poly}(n,L,\delta^{-1}))$, scaling linearly with
the stabilizer extent $\xi_{\mathsf{STAB}_{n}}(\ket{T}^{\otimes t})$ of $t$
copies of the single-qubit magic state $\ket{T}$ associated with a $T$-gate
[11].

###  1.3 The fermionic Gaussian extent

In the following, we generalize the notion of the extent beyond stabilizer
computations to any dictionary $\mathcal{D}$. We refer to the corresponding
quantity as the $\mathcal{D}$-extent $\xi_{\mathcal{D}}(\Psi)$ of $\Psi$. We
assume throughout that we are interested in pure state quantum computations on
a Hilbert space $\mathcal{H}$, and that the dictionary $\mathcal{D}$ is a
subset of pure states on $\mathcal{H}$. Then the $\mathcal{D}$-extent
$\xi_{\mathcal{D}}(\Psi)$ of $\Psi\in\mathcal{H}$ is defined as

$\displaystyle\xi_{\mathcal{D}}(\Psi)$
$\displaystyle=\inf_{N\in\mathbb{N}}\inf_{\varphi_{1},\ldots,\varphi_{N}\in\mathcal{D}}\left\\{\|\gamma\|^{2}_{1}\
\big{|}\ \gamma\in\mathbb{C}^{N}\textrm{ such that
}\sum_{j=1}^{N}\gamma_{j}\varphi_{j}=\Psi\right\\}\ .$ (8)

Here $\|\gamma\|_{1}=\sum_{j=1}^{N}|\gamma_{j}|$ is the $L^{1}$-norm of the
vector $\gamma$. That is, the $\mathcal{D}$-extent $\xi_{\mathcal{D}}(\Psi)$
is the $L^{1}$-norm of the coefficients minimized over all decompositions of
$\Psi$ into a finite linear combination of elements of the dictionary
$\mathcal{D}$. As mentioned above, quantities of the form (8) are well-studied
in the context of signal-processing.

When the dictionary $\mathcal{D}$ is a finite subset of a Hilbert space
$\mathcal{H}\cong\mathbb{C}^{d}$, the $\mathcal{D}$-extent
$\xi_{\mathcal{D}}(\Psi)$ of a state $\Psi\in\mathcal{H}$ can be expressed as
a second-order cone program [28] (see also e.g., [29]), as in Appendix A of
Ref. [19]. Second-order cone programs can be solved in time polynomial in
$\max(d,|\mathcal{D}|)$. We are typically interested in cases where
$\mathcal{D}$ contains a basis of $\mathcal{H}$ (such that every state can
indeed be represented as a linear combination of dictionary elements): Here
this runtime is at least polynomial in $d$. For example, in the case
$\mathcal{D}=\mathsf{STAB}_{n}$ of stabilizer states on $n$ qubits, this leads
to an exponential scaling in $n^{2}$. Beyond algorithmic considerations
related to the evaluation of the extent, the fact that
$\xi_{\mathcal{D}}(\Psi)$ is given by a second-order cone program provides
useful analytical insight by convex programming duality. Indeed, this fact has
previously been exploited both for showing multiplicativity of the stabilizer
extent for states of small dimension [18], as well as to show non-
multiplicativity in high dimensions [19]. In Section 6, we also exploit this
connection to relate the $\mathcal{D}$-fidelity $F_{\mathcal{D}}(\Psi)$ with
the $\mathcal{D}$-extent $\xi_{\mathcal{D}}(\Psi)$.

In contrast, the $\mathcal{D}$-extent $\xi_{\mathcal{D}}(\Psi)$ for an
infinite, i.e., continuously parameterized, dictionary $\mathcal{D}$
constitutes additional mathematical challenges as an optimization problem.
This is the case of interest here as we are considering the dictionary
$\mathcal{D}=\mathcal{G}_{n}$ consisting of all $n$-fermion Gaussian states in
the following. We call the associated quantity $\xi_{\mathcal{G}_{n}}(\Psi)$
the (fermionic) Gaussian extent of an $n$-fermion state $\Psi$. Our focus here
is on discussing the role of the quantity $\xi_{\mathcal{G}_{n}}(\Psi)$ in the
context of classically simulating fermionic linear optics, and its behavior on
tensor products. A detailed discussion of the algorithmic problem of computing
$\xi_{\mathcal{G}_{n}}(\Psi)$ for an arbitrary state $\Psi$, and finding a
corresponding optimal decomposition of $\Psi$ into a linear combination of
Gaussian states is beyond the scope of this work. We refer to e.g., [30] where
semidefinite relaxations are given for the related atomic norm minimization
problem in cases where the atomic set (corresponding to the dictionary) has
algebraic structure. Similar techniques may be applicable to the fermionic
Gaussian extent.

###  1.4 On the (sub)multiplicativity of the extent

Consider a situation where an operation $E\not\in\mathcal{E}$ not belonging to
the set $\mathcal{E}$ of efficiently simulable operations is implemented by
using a “magic” resource state $\Psi\not\in\mathcal{D}$. For example, if
$\mathcal{D}=\mathsf{STAB}_{n}$ is the set of stabilizer states, $\mathcal{E}$
the set of Clifford unitaries and $\mathcal{M}$ the set of single-qubit
Pauli-$Z$-measurements, then a non-Clifford gate (such as the $T$-gate) can be
realized by an (adaptive) Clifford circuit at the cost of consuming a non-
Clifford state (such as the state $\ket{T}$) [11]. Similar “gadget
constructions” exist for fermionic linear optics, where non-Gaussian unitaries
are realized by Gaussian unitaries and non-Gaussian states [12, 13]. A natural
question arising in this situation is to characterize the cost of simulating
the application of two independent magic gates
$E_{1},E_{2}\not\in\mathcal{E}$, each realized by efficiently simulable
operations (belonging to $\mathcal{E}$) using magic states
$\Psi_{1},\Psi_{2}$. For any reasonable simulation algorithm, we expect the
required simulation effort to increase at most multiplicatively. Indeed, this
feature is reflected in the submultiplicativity property

$\displaystyle\xi_{\mathcal{D}_{3}}(\Psi_{1}\otimes\Psi_{2})$
$\displaystyle\leq\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}\ $ (9)

of the $\mathcal{D}$-extent. In Eq. (9), we are considering Hilbert spaces
$\mathcal{H}_{1}$, $\mathcal{H}_{2}$ and their tensor product
$\mathcal{H}_{3}=\mathcal{H}_{1}\otimes\mathcal{H}_{2}$, as well as
dictionaries $\mathcal{D}_{j}\subset\mathcal{H}_{j}$ for $j\in[3]$. The
submultiplicativity property (9) follows immediately from the definition of
the extent if the three dictionaries satisfy the inclusion property

$\displaystyle\mathcal{D}_{1}\otimes\mathcal{D}_{2}\subset\mathcal{D}_{3}\ .$
(10)

In particular, this is satisfied e.g., when the dictionary
$\mathcal{D}_{j}=\mathsf{STAB}_{n_{j}}\subset(\mathbb{C}^{2})^{\otimes n_{j}}$
is the set of $n_{j}$-qubit stabilizer states for $j\in[3]$, with
$n_{3}=n_{1}+n_{2}$, or when considering the set of (even) Gaussian states
(see below).

While the submultiplicativity property (9) is a trivial consequence of Eq.
(10), the question of whether or not the stronger multiplicativity property

$\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}\ $ (11)

holds for the $\mathcal{D}$-extent is a much less trivial problem. If the
multiplicativity property (11) is satisfied, then computing the extent of a
product state can be broken down into several smaller optimization problems:
It suffices to compute the extent of each factor in the tensor product.
Furthermore, the classical simulation cost (with typical algorithms) when
applying several non-free (“magic”) gates constructed by gadgets increases at
an exponential rate determined by the individual gates. In contrast, if the
extent is not multiplicative (i.e., the equality in (11) is not satisfied for
some states $\Psi_{j}\in\mathcal{H}_{j}$, $j\in[2]$), then such a
simplification is not possible. More surprisingly, such a violation of
multiplicativity implies that the classical simulation cost of applying
certain non-free gates can be reduced by treating these jointly instead of
individually. We note that in the slightly different context of so-called
circuit knitting, similar savings in complexity have been shown to be
significant [31].

Previous work established that the stabilizer extent is multiplicative even
for multiple factors, that is,

$\displaystyle\xi_{\mathsf{STAB}_{n_{1}+\cdots+n_{r}}}(\Psi_{1}\otimes\cdots\otimes\Psi_{r})$
$\displaystyle=\prod_{j=1}^{r}\xi_{\mathsf{STAB}_{n_{j}}}(\Psi_{j})\qquad\textrm{
for all }\qquad\Psi_{j}\in(\mathbb{C}^{2})^{\otimes n_{j}},j\in[r]$ (12)

if the factors are single-qubit, $2$\- or $3$-qubit states, i.e.,
$n_{j}\in[3]$, see Ref. [18]. An example is the stabilizer extent of a tensor
product of $t$ copies of the magic (single-qubit) state
$\ket{T}=(\ket{0}+e^{i\pi/4}\ket{1})/\sqrt{2}$ associated with the $T$-gate.
Multiplicativity for qubit states gives
$\xi_{\mathsf{STAB}_{1}}(\ket{T}^{\otimes
t})=\xi_{\mathsf{STAB}_{t}}(\ket{T})^{t}$, where
$\xi_{\mathsf{STAB}_{1}}(\ket{T})$ is known to be approximately $1.17$ [17].
This translates to an overhead exponential in $t$ in the runtime of stabilizer
computations supplemented with $t$ $T$-gates. Surprisingly, the stabilizer
extent has been shown not to be multiplicative (for all pairs of states) in
high dimensions [19].

For (pure) Gaussian states, the Gaussian extent of a $1$-, $2$\- and $3$-mode
pure fermionic state is trivially one because any $1$-, $2$\- and $3$-mode
pure fermionic state is Gaussian [32] and is thus an element of the
dictionary. Hence the Gaussian extent is (trivially) multiplicative if the
factors are $1$-, $2$\- or $3$-mode fermionic states. The simplest non-trivial
case is that of $n=4$ fermionic modes in each factor.

###  1.5 Our contribution

Our results concern fermionic linear optics, the computational model
introduced in Section 1.1.2 described by the triple
$(\mathcal{G}_{n},\mathcal{E}_{\textrm{Gauss}},\mathcal{M}_{\textrm{number}})$
of fermionic Gaussian pure states on $n$ fermions, Gaussian unitary operations
and number state measurements. We propose classical simulation algorithms for
the case where the initial state $\Psi\in\mathcal{H}_{n}$ is an arbitrary pure
state in the $n$-fermion Hilbert space $\mathcal{H}_{n}$ (instead of belonging
to the set $\mathcal{G}_{n}\subset\mathcal{H}_{n}$ of Gaussian states). Our
results are two-fold:

New simulation algorithms.  We give algorithms realizing the functionalities
described in Section 1.1 exactly for the triple
$(\mathcal{H}_{n},\mathcal{E}_{\textrm{Gauss}},\mathcal{M}_{\textrm{number}})$.
This immediately gives rise to efficient algorithms for weak and strong
simulation of circuits with non-Gaussian initial states. The corresponding
runtimes of these building blocks, which we refer to as
$(\chi\mathsf{evolve},\chi\mathsf{measureprob},\allowbreak\chi\mathsf{postmeasure})$,
depend on the Gaussian rank $\chi=\chi_{\mathcal{G}_{n}}(\Psi)$ of the initial
state $\Psi$ and are summarized in Table 3.

algorithm | $\mathsf{time}$
---|---
$\chi\mathsf{evolve}$ | $O(\chi n^{3})$
$\chi\mathsf{measureprob}$ | $O(\chi^{2}n^{3})$
$\chi\mathsf{postmeasure}$ | $O(\chi n^{3})$
Table 3: Runtimes of the building blocks $\chi\mathsf{evolve}$,
$\chi\mathsf{measureprob}$, $\chi\mathsf{postmeasure}$ for exact simulation of
$n$-qubit fermionic linear optics circuits with a non-Gaussian initial state
$\Psi$ of Gaussian rank $\chi=\chi_{\mathcal{G}_{n}}(\Psi)$. Evolution
corresponds to the application of a Gaussian unitary from a set
$\mathcal{E}_{\textrm{Gauss}}$ of generators (specified below), and the set of
measurements is given by occupation number measurements on each of the modes.

Key to the construction of these algorithms is a novel way of keeping track of
relative phases in superpositions of Gaussian states, see Section 3. We argue
that our techniques can be applied more generally to adapt simulation
procedures developed, e.g., for Clifford circuits, to the setting of fermionic
linear optics. In order to illustrate this procedure, we apply it to the
simulation algorithms of [17, 18] for efficient (approximate) classical
simulation algorithms. In this way, we obtain new approximate simulation
algorithms with runtimes depending linearly on the fermionic Gaussian extent
$\xi=\xi_{\mathcal{G}_{n}}(\Psi)$ of the initial state $\Psi$, see Table 4 for
a summary of the corresponding runtimes.

algorithm | $\mathsf{time}$
---|---
$\mathsf{approx}\mathsf{evolve}$ | $O(\xi\delta^{-2}n^{3})$
$\mathsf{approx}\mathsf{measureprob}$ | $O(\xi\delta^{-2}\epsilon^{-2}p_{f}^{-1}n^{7/2})$
$\mathsf{approx}\mathsf{postmeasure}$ | $O(\xi\delta^{-2}n^{3})$
Table 4: Runtimes of building blocks $\mathsf{approx}\mathsf{evolve}$,
$\mathsf{approx}\mathsf{measureprob}$, $\mathsf{approx}\mathsf{postmeasure}$
for approximate simulation of $n$-qubit fermionic linear optics circuits with
a non-Gaussian initial state $\Psi$ of Gaussian extent
$\xi=\xi_{\mathcal{G}_{n}}(\Psi)$. The parameters $(\epsilon,\delta,p_{f})$
determine the quality of the approximation.

They depend inverse-polynomially on parameters $(\delta,\epsilon,p_{f})$
determining the accuracy of the simulation. The error $\delta$ describes a
certain “offset”, i.e., a systematic error: Instead of simulating the dynamics
of the circuit with the (ideal) initial state $\Psi$, the simulation algorithm
emulates the dynamics when using a different starting state $\tilde{\Psi}$
which is $\delta$-close to $\Psi$, i.e., which satisfies
$\|\Psi-\tilde{\Psi}\|\leq\delta$. The algorithm
$\mathsf{approx}\mathsf{evolve}$ computes evolution exactly on the state used
in the simulation (i.e., it preserves the approximation error $\delta$
relative to the ideal initial state). In contrast, the procedure
$\mathsf{approx}\mathsf{measureprob}$ can fail with probability $p_{f}$, and
both $\mathsf{approx}\mathsf{measureprob}$ and
$\mathsf{approx}\mathsf{postmeasure}$ introduce an additional error quantified
by $\epsilon$ (if $\mathsf{approx}\mathsf{measureprob}$ succeeds): Instead of
returning the probability $p(0)$ of obtaining zero occupation number when
measuring the state, the output of $\mathsf{approx}\mathsf{measureprob}$ is a
value $\tilde{p}$ which satisfies $|\tilde{p}-p(0)|\leq O(\epsilon)$.
Similarly, the output of $\mathsf{approx}\mathsf{postmeasure}$ is a
description of a state that is $O(\epsilon)$-close to the actual post-
measurement state. These parameters and runtimes are analogous to those
obtained in [18] for simulating Clifford circuits with non-stabilizer initial
states. In particular, they imply that a circuit with initial state $\Psi$
involving $T$ Gaussian unitaries and $L$ occupation number measurements can be
weakly simulated in time $\tilde{O}(\epsilon^{-2}\xi)$, such that the sampled
measurement outcomes are $\epsilon$-close in $L^{1}$-distance to the ideal
(joint) output distribution of all measurements. Here the notation
$\tilde{O}(\cdot)$ suppresses a factor polynomial in $n,T,L$ and
$\log(\epsilon^{-1})$, see [17] for details.

On the multiplicativity of the Gaussian extent and the Gaussian fidelity.
Motivated by the relevance of the Gaussian extent
$\xi_{\mathcal{G}_{n}}(\Psi)$ for characterizing the complexity of classical
simulation, we study multiplicativity properties of both the
$\mathcal{D}$-fidelity $F_{\mathcal{D}}(\Psi)$ as we well as the
$\mathcal{D}$-extent $\xi_{\mathcal{D}}(\Psi)$ for a general infinite, i.e.,
continuously parameterized, dictionary $\mathcal{D}$. We show that
multiplicativity of the $\mathcal{D}$-fidelity is closely related to that of
the $\mathcal{D}$-extent: For a general family of (discrete or continuous)
dictionaries $\mathcal{D}_{j}\subset\mathcal{H}_{j}$ for $j\in[3]$ with the
property

$\displaystyle\mathcal{D}_{1}\otimes\mathcal{D}_{2}\subset\mathcal{D}_{3}\ ,$
(13)

multiplicativity of the $\mathcal{D}$-fidelity, i.e.,

$\displaystyle
F_{\mathcal{D}_{3}}\left(\Psi_{1}\otimes\Psi_{2}\right)=F_{\mathcal{D}_{1}}\left(\Psi_{1}\right)F_{\mathcal{D}_{2}}\left(\Psi_{2}\right)\quad\text{
for all }\quad\Psi_{j}\in\mathcal{H}_{j}\text{ for }j\in[2]$ (14)

implies multiplicativity of the $\mathcal{D}$-extent, i.e.

$\displaystyle\xi_{\mathcal{D}_{3}}\left(\Psi_{1}\otimes\Psi_{2}\right)=\xi_{\mathcal{D}_{1}}\left(\Psi_{1}\right)\xi_{\mathcal{D}_{2}}\left(\Psi_{2}\right)\quad\text{
for all }\quad\Psi_{j}\in\mathcal{H}_{j}\text{ for }j\in[2]\ .$ (15)

We note that for stabilizer states $\mathcal{D}=\mathsf{STAB}_{n}$, a similar
route was followed in Ref. [18] to show multiplicativity of the stabilizer
extent $\xi_{\mathsf{STAB}_{n}}$ with respect to the tensor product of $1$-,
$2$\- and $3$-qubit states. Our main contribution is an extension of this
connection to the case of infinite dictionaries by the use of nets. We expect
this connection to be helpful in proving or disproving multiplicativity of the
extent more generally.

We subsequently make use of this connection to the Gaussian fidelity to show
that the fermionic Gaussian extent is multiplicative for the tensor product of
any two 4-mode fermionic states with positive parity, i.e.,

$\displaystyle\xi_{\mathcal{G}_{8}}(\Psi_{1}\otimes\Psi_{2})=\xi_{\mathcal{G}_{4}}(\Psi_{1})\xi_{\mathcal{G}_{4}}(\Psi_{2})\qquad\textrm{
for all }\qquad\Psi_{1},\Psi_{2}\in\mathcal{H}^{4}_{+}\ .$ (16)

Here $\mathcal{H}^{4}_{+}$ denotes the set of $4$-mode fermionic states with
positive parity. The proof of (16) relies on the Schmidt decomposition of
fermionic Gaussian states and specific properties of $4$-mode (positive
parity) fermionic states.

The result (16) gives the first non-trivial example of multiplicativity of the
Gaussian extent. Multiplicativity for more general cases such as that of
multiple $4$-mode fermionic factors remains an open problem.

###  1.6 Prior and related work

The starting point of our work is the fact that fermionic Gaussian operations
acting on Gaussian states can be efficiently simulated classically as shown in
pioneering work by Terhal and DiVincenzo [9] and Knill [10]. The model and its
simulability are closely related to that of matchgate computations introduced
by Valiant [8], where so-called matchgates correspond to a certain certain
subset of Gaussian operations (see also [33]). In analogy to the fermionic
context, the efficient simulability of bosonic Gaussian circuits was
recognized at around the same time [34, 35]. In an effort to identify
commonalities between simulation algorithms for a variety of quantum
computational models, Somma et al. [36] provided a unifying Lie algebraic
treatment which gives a counterpart to the Gottesman-Knill theorem for the
simulability of Clifford circuits [37, 7].

While matchgate circuits, fermionic and bosonic linear optics, and Clifford
circuits provide rich classes of efficiently simulable models for the study of
many-body dynamics associated with quantum circuits, it is desirable to extend
the applicability of such simulation methods. There has been significant
interest in this problem resulting in a range of approaches. We only briefly
discuss these here to give an overview, without attempting to give an
exhaustive treatment.

A first prominent approach is the use of quasi-probability distributions to
describe states and corresponding dynamics. Such a description typically
applies to a subset of density operators: For example, it has been shown in
[38, 39] in the context of continuous-variable systems that circuits applying
bosonic Gaussian operations to initial states with a positive Wigner function
(a strict superset of the set of bosonic Gaussian states) can be simulated
efficiently. Negativity of the Wigner function (both in the continuous-
variable as well as the qubit context) thus serves as a resource for quantum
computation, also see e.g., [40, 41]. It is also closely related to
contextuality, see [42], and thus connects contextuality to the complexity of
classical simulation [43, 44]. Not unlike the notorious sign problem in
quantum Monte-Carlo methods applied in many-body physics, the runtimes of
corresponding (randomized) simulation algorithms scale with certain measures
of “negativity” of the initial state.

The concept of a convex-Gaussian state was introduced and studied in [32] to
extend the range of fermionic linear optics simulation methods. This is
related to quasi-probability representations in the sense that initial states
of a particular form are shown to lead to efficient simulability. Here a
density operator is called convex-Gaussian if it is a convex combination of
fermionic Gaussian states. The utility of this concept was illustrated in [32]
by showing a converse to the fault-tolerance threshold theorem: Sufficiently
noisy quantum circuits can be simulated classically because the corresponding
states turn out to be convex-Gaussian. A detailed characteriziation of convex-
Gaussianity is necessary to translate this into explicit (numerical) threshold
estimates. An infinite hierarchy of semidefinite programs was constructed in
[32] to detect convex-Gaussianity, and this was subsequently shown to be
complete [45]. This hierarchy also provides a way of determining whether a
state is close to being convex-Gaussian [45].

A second important class of approaches are rank-based methods. Here the non-
free resource (either a state or an operation) is decomposed into a linear
combination of free (i.e., efficiently simulable) resources. Our work follows
this approach, which is detailed in Section 4.1 for pure states. For Clifford
computations, this involves writing general states as superpositions of
stabilizer states. The development of such simulators was pioneered by Bravyi,
Smith, and Smolin [16] with subsequent work dealing with approximate
stabilizer decompositions [17].

The concept of low-rank (approximate) decompositions of quantum states or
operations into more easily treatable basic objects appears in a variety of
forms: For example, the work [18] also discusses – in addition to state vector
decompositions – decompositions of non-Clifford unitaries into sums of
Clifford operations. In Ref. [46], a similar approach was taken to
approximately decompose non-Gaussian fermionic unitary operations into linear
combinations of Gaussian channels. In all these cases, the main challenge is
to identify optimal (or simply good) decompositions (e.g., in terms of rank or
an extent-like quantity).

In more recent work, Mocherla, Lao and Browne [47] study the problem of
simulating matchgate circuits using universality-enabling gates. They provide
a simulation algorithm and associated runtime estimates for estimating
expectation values of single-qubit observables in output states obtained by
applying a matchgate circuit to a product state. This problem is closely
related to the problem considerd in this work as matchgate circuits
efficiently describe evolution under a quadratic fermionic Hamiltonian. The
approach taken in [47] is quite different from ours, however: The classical
simulator keeps track of the density operator by tracking its coefficients in
the Pauli (operator) basis, using the structure of corresponding linear maps
associated with matchgates. The effect of a specific set of universality-
enabling gates is then analyzed in detail. This extends the sparse simulation
method for matchgate circuits to circuits augmented with such gates. The
runtime estimates of [47] apply to certain universality-providing gates. In
contrast, our constructions can in principle also be applied to (gadget-based)
constructions of arbitrary gates and provide gate-specific information. For
gates close to the identity, for example, this may provide additional resource
savings (in terms of e.g., the rate of growth for several uses of such a
gate).

Near the completion of our work, we became aware of concurrent independent
work on fermionic circuits with non-Gaussian operations, see the papers [48,
49] which were posted simultaneously with our work to the arXiv. Reference
[48] proposes an alternative classical simulation algorithm for such circuits
whose basic building blocks have identical runtime as our algorithm. In
particular, when applying a generator $\exp(\vartheta/2c_{j}c_{k})$ of the
Gaussian unitary group, the runtime of the corresponding state update routine
is $O(n^{3})$ as in our work. Specifying a circuit in terms of a sequence of
generators is a common assumption when discussing restricted gate sets. If
instead, a general Gaussian unitary is specified by an element $R\in SO(2n)$,
the runtimes are different: By decomposing $R$ into $O(n^{2})$ Givens
rotations (generators), our techniques yield a runtime of order $O(n^{5})$,
whereas the algorithm of [48] runs in time $O(n^{4})$ in this case. A more
significant distinction is the underlying approach: While we rely on a minimal
extension of the covariance matrix formalism, the approach of [48] involves an
explicit decomposition of fermionic Gaussian unitaries inspired by a certain
canonical form of Clifford circuits [18]. This decomposition of a Gaussian
unitary $U$ is subsequently used to specify a Gaussian state $U|0\rangle$
including its global phase.

The work [49] analyses the set of fermionic Gaussian states to investigate
properties of the (exact and approximate) Gaussian rank, the Gaussian fidelity
and the Gaussian extent. A major result obtained in [49] is the
multiplicativity of the extent for a tensor product state with any number of
$4$-mode factors. This mirrors what is known for two-qubit and three-qubit
factors in the stabilizer context [18]. The proof in Ref. [49] is inspired by
Ref. [19] and makes use of tools from convex optimization. It involves the use
of a net of Gaussian states to establish an upper bound on the Gaussian
fidelity. While we only establish multiplicativity for two factors (i.e., a
special case of their result), our approach is different and may be of
independent interest: We use a net to show that multiplicativity of the
Gaussian fidelity implies that of the Gaussian extent, and employ the Schmidt
decomposition [50] of fermionic Gaussian states.

###  1.7 Outline

The paper is structured as follows. In Section 2, we give some background on
fermionic linear optics, reviewing fermionic Gaussian operations and states,
inner product formulas for Gaussian states and tensor products of fermionic
systems. In Sections 3 and 4 we describe classical algorithms for simulation
of Gaussian and non-Gaussian fermionic circuits, respectively. Specifically,
in Section 3 we provide an algorithm $\mathsf{overlap}$ for computing the
overlap of two Gaussian states, an algorithm $\mathsf{evolve}$ to simulate
unitary evolution of a Gaussian state, and algorithms $\mathsf{measureprob}$
and $\mathsf{postmeasure}$ to simulate measurements of occupation numbers. All
these algorithms keep track of the phase of the state. In Section 4 we extend
the simulation described in Section 3 to allow for non-Gaussian input states.
The remainder of this work is focused on the multiplicativity of the fermionic
Gaussian extent. In Section 5, we prove the multiplicativity of the fermionic
Gaussian fidelity for the tensor product of any two 4-mode fermionic states
with positive parity. Section 6 is devoted to showing that the
multiplicativity of the $\mathcal{D}$-fidelity implies multiplicativity of the
$\mathcal{D}$-extent for general (finite and infinite, i.e., continuously
parameterized) dictionaries. Finally, the results from Sections 5 and 6 are
used to prove the main result in Section 7, namely the multiplicativity of the
fermionic Gaussian extent for the tensor product of any two 4-mode fermionic
states with positive parity.

##  2 Background

In this section, we give some background on fermionic linear optics to fix
notation.

###  2.1 Dirac and Majorana operators

Throughout, we consider fermionic systems composed of $n$ modes, with (Dirac)
creation- and annihilation operators $a_{j}^{\dagger},a_{j}$, $j\in[n]$,
satisfying the canonical anticommutation relations

$\displaystyle\\{a_{j},a_{k}^{\dagger}\\}=\delta_{j,k}I\qquad\textrm{ and
}\qquad\\{a_{j},a_{k}\\}=\\{a_{j}^{\dagger},a_{k}^{\dagger}\\}=0\qquad\textrm{
for all }\qquad j,k\in[n]\ .$ (17)

The fermionic vacuum state $\ket{0_{F}}$ is the (up to a phase) unique unit
vector satisfying $a_{j}\ket{0_{F}}=0$ for all $j\in[n]$. For
$x=(x_{1},\ldots,x_{n})\in\\{0,1\\}^{n}$, we define the number state $\ket{x}$
by

$\displaystyle\ket{x}=(a_{1}^{\dagger})^{x_{1}}\cdots(a_{n}^{\dagger})^{x_{n}}\ket{0_{F}}\
.$ (18)

The states $\\{\ket{x}\\}_{x\in\\{0,1\\}^{n}}$ are an orthonormal basis of the
underlying Hilbert space $\mathcal{H}^{n}\cong(\mathbb{C}^{2})^{\otimes n}$. A
state $\ket{x}$ is a simultaneous eigenstate of the occupation number
operators $a_{j}^{\dagger}a_{j}$, $j\in[n]$, where $x_{j}$ is the eigenvalue
of $a_{j}^{\dagger}a_{j}$. For later reference, we note that

$\displaystyle a_{j}\ket{x}=(-1)^{\eta_{j}(x)}x_{j}\ket{x\oplus
e_{j}}\qquad\text{and}\qquad
a^{\dagger}_{j}\ket{x}=(-1)^{\eta_{j}(x)}\overline{x_{j}}\ket{x\oplus e_{j}}\
,$ (19)

with the definition

$\displaystyle\eta_{j}(x)=\sum_{k=1}^{j-1}x_{k}\qquad\textrm{ for }\qquad
j\in[n]\ ,$ (20)

where we write $\overline{0}=1$ and $\overline{1}=0$, where
$e_{j}\in\\{0,1\\}^{n}$ is given by $(e_{j})_{k}=\delta_{j,k}$ for $k\in[n]$,
and where $\oplus$ denotes bitwise addition modulo $2$.

It will be convenient to work with Majorana operators $\\{c_{j}\\}_{j=1}^{2n}$
defined by

$\displaystyle c_{2j-1}=a_{j}+a_{j}^{\dagger}\qquad\text{and}\qquad
c_{2j}=i\left(a_{j}-a_{j}^{\dagger}\right)\ .$ (21)

Majorana operators are self-adjoint and satisfy the relations

$\displaystyle\\{c_{j},c_{k}\\}=2\delta_{jk}I\qquad\textrm{ and }\qquad
c_{j}^{2}=I\qquad\textrm{ for }\qquad j,k\in[2n]\ .$ (22)

For $\alpha\in\\{0,1\\}^{2n}$, we call the self-adjoint operator

$\displaystyle c(\alpha)$
$\displaystyle=i^{|\alpha|\cdot(|\alpha|-1)/2}c_{1}^{\alpha_{1}}\cdots
c_{2n}^{\alpha_{2n}}\ $ (23)

a Majorana monomial. Here $|\alpha|=\sum_{j=1}^{2n}\alpha_{j}$ denotes the
Hamming weight of $\alpha\in\\{0,1\\}^{2n}$. The set
$\\{c(\alpha)\\}_{\alpha\in\\{0,1\\}^{n}}$ constitutes an orthonormal basis of
the real vector space of self-adjoint operators on $\mathcal{H}^{n}$ equipped
with the (normalized) Hilbert-Schmidt inner product $\langle
A,B\rangle=2^{-n}\mathop{\mathrm{tr}}\nolimits(A^{\dagger}B)$. The Majorana
monomials satisfy

$\displaystyle c(y)c(x)=(-1)^{|x|\cdot|y|+x\cdot
y}c(x)c(y)\qquad\text{with}\qquad x,y\in\\{0,1\\}^{2n}\ ,$ (24)

where $x\cdot y=\sum_{j=1}^{n}x_{j}y_{j}$. In particular, if either $x$ or $y$
have even Hamming weight then $c(x)c(y)=(-1)^{x\cdot y}c(y)c(x)$. In the
following, we will denote the set of even- and odd-weight $2n$-bit strings by
$\\{0,1\\}_{+}^{2n}$ and $\\{0,1\\}_{-}^{2n}$, respectively.

The parity operator

$\displaystyle P=i^{n}c_{1}c_{2}\cdots c_{2n}\ $ (25)

is the Majorana monomial associated with $\alpha=1^{2n}=(1,\ldots,1)$. The
parity operator commutes with every even-weight Majorana monomial and anti-
commutes with every odd-weight Majorana monomial, i.e., we have

$\displaystyle Pc(\alpha)$
$\displaystyle=(-1)^{|\alpha|}c(\alpha)P\qquad\textrm{ for every
}\qquad\alpha\in\\{0,1\\}^{2n}\ .$ (26)

The Hilbert space
$\mathcal{H}^{n}=\mathcal{H}^{n}_{+}\oplus\mathcal{H}^{n}_{-}$ associated with
$n$ fermions decomposes into a direct sum of positive- and negative-parity
vectors

$\displaystyle\mathcal{H}^{n}_{+}$
$\displaystyle=\left\\{\Psi\in\mathcal{H}^{n}\ |\ P\Psi=\Psi\right\\}\ \text{
and }$ (27) $\displaystyle\mathcal{H}^{n}_{-}$
$\displaystyle=\left\\{\Psi\in\mathcal{H}^{n}\ |\ P\Psi=-\Psi\right\\}\ .$
(28)

We call a state $\Psi\in\mathcal{H}^{n}$ of definite parity if either
$\Psi\in\mathcal{H}^{n}_{+}$ or $\Psi\in\mathcal{H}^{n}_{-}$. An element
$X\in\mathcal{B}(\mathcal{H}^{n})$ belonging to the set
$\mathcal{B}(\mathcal{H}^{n})$ of linear operators on $\mathcal{H}^{n}$ is
called even (odd) if it is a linear combination of Majorana monomials
$c(\alpha)$ with $\alpha\in\\{0,1\\}^{2n}$ of even (odd) weight. An immediate
consequence of these definitions is that a state $\Psi\in\mathcal{H}^{n}$ has
definite parity if and only if $|\Psi\rangle\langle\Psi|$ is even (see e.g.,
[51, Proposition 1] for a proof).

###  2.2 Gaussian unitaries

A unitary operator $U$ on $\mathcal{H}^{n}$ is Gaussian if and only if it maps
a Majorana operator $c_{j}$ to a linear combination of Majorana operators,
i.e.

$\displaystyle Uc_{j}U^{\dagger}=\sum_{k=1}^{2n}R_{jk}c_{k}\ ,$ (29)

where $R\in O(2n)$ is a real orthogonal matrix. Ignoring overall phases, the
group of Gaussian unitary operators is generated by operators of the form

$\displaystyle U_{j,k}(\vartheta)=\exp(\vartheta/2c_{j}c_{k})\qquad\textrm{
with }\qquad\vartheta\in[0,2\pi)\textrm{ and }j<k\in[2n]$ (30)

and by operators

$\displaystyle U_{j}=c_{j}\qquad\textrm{ with }\qquad j\in[2n]\ .$ (31)

The operator $U_{j,k}(\vartheta)$ implements the rotation

$\displaystyle U_{j,k}(\vartheta)c_{j}U_{j,k}(\vartheta)^{\dagger}$
$\displaystyle=\cos(\vartheta)c_{j}-\sin(\vartheta)c_{k}$ (32) $\displaystyle
U_{j,k}(\vartheta)c_{k}U_{j,k}(\vartheta)^{\dagger}$
$\displaystyle=\sin(\vartheta)c_{j}+\cos(\vartheta)c_{k}$ $\displaystyle
U_{j,k}(\vartheta)c_{\ell}U_{j,k}(\vartheta)^{\dagger}$
$\displaystyle=c_{\ell}$ $\displaystyle\qquad\textrm{for
}\qquad\ell\not\in\\{j,k\\}\ .$

The operator $U_{j}=c_{j}$ leaves $c_{j}$ invariant and flips the sign of each
$c_{k}$ with $k\neq j$, i.e., it implements the reflection

$\displaystyle U_{j}c_{j}U_{j}^{\dagger}$ $\displaystyle=c_{j}$ (33)
$\displaystyle U_{j}c_{k}U_{j}^{\dagger}$ $\displaystyle=-c_{k}\qquad\textrm{
for }\qquad k\neq j\ .$ (34)

We note that by relation (26), every generator $U_{j,k}(\vartheta)$ is parity-
preserving, whereas every generator $U_{j}$ reverses the parity, i.e.,

$\displaystyle U_{j,k}(\vartheta)PU_{j,k}(\vartheta)^{\dagger}$
$\displaystyle=P$ $\displaystyle\qquad\text{ for all }\qquad
k>j\in[n],\vartheta\in[0,2\pi)\ ,$ (35) $\displaystyle U_{j}PU_{j}^{\dagger}$
$\displaystyle=-P$ $\displaystyle\qquad\text{ for all }\qquad j\in[n]\ .$

Every orthogonal matrix $R$ gives rise to a Gaussian unitary $U_{R}$
satisfying (29). The unitary $U_{R}$ is unique up to a global phase, and
$R\mapsto U_{R}$ is called the metaplectic representation. We can fix the
global phase of $U_{R}$ uniquely, e.g., by the following procedure. Every
element $R\in O(2n)$ can be uniquely decomposed into a product

$\displaystyle R$ $\displaystyle=S_{0}S_{1}\cdots S_{L}$ (36)

with $L\leq\frac{2n(2n-1)}{2}$ and where

$\displaystyle S_{0}$ $\displaystyle=\begin{cases}I\qquad\textrm{ if }R\in
SO(2n)\\\ R_{1}\qquad\textrm{ otherwise }\ ,\end{cases}$ (37)

where for each $r\in[L]$, the matrix $S_{r}$ is of the form

$\displaystyle S_{r}$
$\displaystyle=R_{j_{r},k_{r}}(\vartheta_{r})\qquad\textrm{ for some }\qquad
j_{r}<k_{r}\in[2n],\vartheta_{r}\in[0,2\pi)\ .$ (38)

Here $R_{1}\in O(2n)$ is associated with the unitary $U_{1}$ by Eq. (34),
whereas $R_{j,k}(\vartheta)\in SO(2n)$ is associated with $U_{j,k}(\vartheta)$
according to Eq. (32). We note that $R_{j,k}(\vartheta)\in SO(2n)$ is a so-
called Givens rotation, introduced in Ref. [52], and a decomposition of the
form (36) can be found by a deterministic algorithm with runtime $O(n^{3})$
(see e.g., Section 5.2.3 in Ref. [53]). In particular, application of this
algorithm defines a function taking an arbitrary element $R\in O(2n)$ to a
unique product of the form (36). Given the (unique) decomposition (36) of
$R\in O(2n)$, we can then define $U_{R}$ as the product

$\displaystyle U_{R}$ $\displaystyle=U_{1}U_{j_{1},k_{1}}(\vartheta_{1})\cdots
U_{j_{L},k_{L}}(\vartheta_{L})\ .$ (39)

Overall, this defines a function $R\mapsto U_{R}$ from $O(2n)$ to the set of
Gaussian unitaries, fixing the phase ambiguity. Throughout the remainder of
this work, $U_{R}$ will denote the Gaussian unitary uniquely fixed by $R$.

###  2.3 Fermionic Gaussian (pure) states

The set of pure fermionic Gaussian states is the orbit of the vacuum state
$\ket{0_{F}}$ under the action of $O(2n)$ defined by the metaplectic
representation, i.e., fermionic Gaussian states are of the form
$U_{R}\ket{0_{F}}$ with $U_{R}$ a fermionic Gaussian unitary. In more detail,
every fermionic Gaussian state $e^{i\theta}U_{R}\ket{0_{F}}$ is uniquely
specified by a pair $(\theta,R)$ with $\theta\in[0,2\pi)$ and $R\in O(2n)$. We
will denote the set of all fermionic Gaussian states by

$\displaystyle\mathcal{G}_{n}=\left\\{e^{i\theta}U_{R}\ket{0_{F}}\ |\
\theta\in[0,2\pi),R\in O(2n)\right\\}\ .$ (40)

By Eq. (35) and because $P\ket{0_{F}}=\ket{0_{F}}$, every pure fermionic
Gaussian state $\Psi$ has a fixed parity, i.e., it is an eigenvector of the
parity operator $P$. This defines a disjoint partition
$\mathcal{G}_{n}=\mathcal{G}_{n}^{+}\cup\mathcal{G}_{n}^{-}$ of the set of
fermionic Gaussian states into positive- and negative-parity states.

###  2.4 Gaussianity condition

In Ref. [54] Bravyi established a necessary and sufficient condition to
determine if a (possibly mixed) state $\rho\in\mathcal{B}(\mathcal{H}^{n})$ is
Gaussian (see Theorem 1 therein). Here Gaussianity of a density operator
$\rho$ is defined by the condition that $\rho$ has the form

$\displaystyle\rho$
$\displaystyle=K\exp\left(i\sum_{j,k=1}^{2n}A_{j,k}c_{j}c_{k}\right)$ (41)

for an antisymmetric matrix $A=-A^{T}\in\mathsf{Mat}_{2n\times
2n}(\mathbb{R})$ and a constant $K>0$. We note that a pure state
$\Psi\in\mathcal{H}^{n}$ is Gaussian if and only if the associated density
operator $\rho=|\Psi\rangle\langle\Psi|$ is Gaussian. (This follows from the
fact that $|0_{F}\rangle\langle
0_{F}|=\frac{1}{2^{n/2}}\exp\left(i\frac{\pi}{4}\sum_{j=1}^{n}c_{2j-1}c_{2j}\right)$.
Indeed, if $\rho=|\Psi\rangle\langle\Psi|$ is a rank-one projection of the
form (41), then it follows from Williamson’s normal form for antisymmetric
matrices that there is $R\in O(2n)$ such that
$RAR^{T}=\bigoplus_{j=1}^{n}\left(\begin{smallmatrix}0&1\\\
-1&0\end{smallmatrix}\right)$. This implies that
$U_{R}|\Psi\rangle\langle\Psi|U_{R}^{\dagger}=|0_{F}\rangle\langle 0_{F}|$ and
thus $\ket{\Psi}=U^{\dagger}_{R}\ket{0_{F}}$. Conversely, for a Gaussian state
$\ket{\Psi}=e^{i\theta}U_{R}\ket{0_{F}}\in\mathcal{G}_{n}$, we can use the
expression for $|0_{F}\rangle\langle 0_{F}|$ to argue that
$|\Psi\rangle\langle\Psi|$ is of the form (41), i.e., Gaussian.) The
characterization of Gaussian density operators established in [54] is the
following.

######  Theorem 2.1 (Theorem 1 in [54]).

Define $\Lambda=\sum_{j=1}^{2n}c_{j}\otimes c_{j}$. An even state
$\rho\in\mathcal{B}(\mathcal{H}^{n})$ is Gaussian if and only if
$[\Lambda,\rho\otimes\rho]=0$.

Based on this characterization [54], the following was shown in [32].

######  Lemma 2.2 (Corollary 1 in [32]).

Let $\rho\in\mathcal{B}(\mathcal{H}^{n})$ be an even state. Then $\rho$ is a
Gaussian pure state if and only if $\Lambda\left(\rho\otimes\rho\right)=0$.

In the following, we only use the statement of Lemma 2.2 applied to pure
states in order to distinguish between Gaussian and non-Gaussian pure states.
We formulate this as follows:

######  Lemma 2.3.

Let $\Psi\in\mathcal{H}^{n}$ be a pure state with fixed parity. Then $\Psi$ is
Gaussian if and only if

$\displaystyle\Lambda(\ket{\Psi}\otimes\ket{\Psi})$ $\displaystyle=0\ .$ (42)

###### Proof.

This follows immediately from the equivalence of the concepts of Gaussianity
of pure states (vectors) and density operators because the density operator
$|\Psi\rangle\langle\Psi|$ is even for any fixed-parity state $\Psi$. ∎

We note that there is an elegant representation-theoretic interpretation of
this characterization of Gaussianity [55]. It is derived from the fact that
Gaussian states are the orbit of the vacuum state $\ket{0_{F}}$ (a highest
weight state) under the action of the metaplectic group, cf. [56, Section IV]
and [57]. We use a version of this reformulation for $4$ fermions, see Lemma
5.1 below, that has first been obtained in [58].

###  2.5 Covariance matrices, Gaussian states and Wick’s theorem

The covariance matrix $\Gamma=\Gamma(\Psi)\in\mathsf{Mat}_{2n\times
2n}(\mathbb{R})$ of a state $\Psi\in\mathcal{H}^{n}$ is the antisymmetric
matrix with entries

$\displaystyle\Gamma_{j,k}(\Psi)$
$\displaystyle=\begin{cases}\langle\Psi,ic_{j}c_{k}\Psi\rangle&\quad\text{for
}j\neq k\\\ 0&\quad\text{for }j=k\end{cases}$ (43)

with $j,k\in[2n]$. It satisfies $\Gamma\Gamma^{T}=I$ for any state
$\Psi\in\mathcal{H}^{n}$.

The expectation value of a Hermitian operator with respect to a Gaussian state
$\Psi$ is fully determined by its covariance matrix $\Gamma=\Gamma(\Psi)$.
This is because the expectation value of a Majorana monomial $c(\alpha)$,
$\alpha\in\\{0,1\\}^{2n}$, is given by Wick’s theorem

$\displaystyle\langle\Psi,c(\alpha)\Psi\rangle=\begin{cases}\mathsf{Pf}(\Gamma[\alpha])\qquad&\textrm{
if }|\alpha|\textrm{ is even }\\\ 0&\textrm{ otherwise }\end{cases}\ .$ (44)

Here $\Gamma[\alpha]\in\mathsf{Mat}_{|\alpha|\times|\alpha|}(\mathbb{R})$ is
the submatrix of $\Gamma$ which includes all rows and columns with index
$j\in[2n]$ such that $\alpha_{j}=1$. Evaluating such expectation values, i.e.,
computing Pfaffians of $|\alpha|\times|\alpha|$-matrices (with $|\alpha|$
even), takes time $O(|\alpha|^{3})$. (Here and below we use the number of
elementary arithmetic operations to quantify the time complexity of
algorithms.)

###  2.6 Inner product formulas for Gaussian states

The modulus of the inner product of two Gaussian states $\Phi_{1},\Phi_{2}$
with identical parity $\sigma\in\\{\pm 1\\}$ and covariance matrices
$\Gamma_{1},\Gamma_{2}$ is given by the expression [59]

$\displaystyle|\langle\Phi_{1},\Phi_{2}\rangle|^{2}$ $\displaystyle=\sigma
2^{-n}\mathsf{Pf}(\Gamma_{1}+\Gamma_{2})\ .$ (45)

For three Gaussian states $\Phi_{0},\Phi_{1},\Phi_{2}$, the expression
$\left\langle\Phi_{0},\Phi_{1}\right\rangle\cdot\left\langle\Phi_{1},\Phi_{2}\right\rangle\cdot\langle\Phi_{2},\Phi_{0}\rangle$
is invariant under a change of the global phase of any of the states, and can
therefore be computed by the covariance matrix formalism. An explicit
expression was derived by Löwdin in [59]. In Ref. [60] Bravyi and Gosset gave
the formula

$\displaystyle\left\langle\Phi_{0},\Phi_{1}\right\rangle\cdot\left\langle\Phi_{1},\Phi_{2}\right\rangle\cdot\left\langle\Phi_{2},\Phi_{0}\right\rangle=\sigma
4^{-n}i^{n}\operatorname{Pf}\left(\begin{array}[]{ccc}i\Gamma_{0}&-I&I\\\
I&i\Gamma_{1}&-I\\\ -I&I&i\Gamma_{2}\end{array}\right)$ (49)

for three Gaussian states $\\{\Phi_{j}\\}_{j=0}^{2}$ of identical parity
$\sigma\in\\{\pm 1\\}$, where $\Gamma_{j}=\Gamma(\Phi_{j})$ is the covariance
matrix of $\Phi_{j}$ for $j=0,1,2$. More generally, they obtained the formula

$\displaystyle\left\langle\Phi_{0},\Phi_{1}\right\rangle\cdot\left\langle\Phi_{1},c(\alpha)\Phi_{2}\right\rangle\cdot\left\langle\Phi_{2},\Phi_{0}\right\rangle=\sigma
4^{-n}i^{n+|\alpha|\cdot(|\alpha|-1)/2}\operatorname{Pf}\left(R_{\alpha}\right)\
,$ (50)

for any even-weight Majorana monomial $c(\alpha)$,
$\alpha\in\\{0,1\\}^{2n}_{+}$, where

$\displaystyle R_{\alpha}=\left(\begin{array}[]{cccc}i\Gamma_{0}&-I&I&0\\\
I&i\Gamma_{1}&-I&0\\\
-I&I&iD_{\alpha}\Gamma_{2}D_{\alpha}&J_{\alpha}^{T}+iD_{\alpha}\Gamma_{2}J_{\alpha}^{T}\\\
0&0&-J_{\alpha}+iJ_{\alpha}\Gamma_{2}D_{\alpha}&iJ_{\alpha}\Gamma_{2}J_{\alpha}^{T}\end{array}\right)\in\mathsf{Mat}_{(6n+|\alpha|)\times(6n+|\alpha|)}(\mathbb{R})\
.$ (55)

Here $D_{\alpha}=\mathsf{diag}(\\{1-\alpha_{j}\\}_{j=1}^{2n})$ is a diagonal
matrix, whereas $J_{\alpha}\in\mathsf{Mat}_{|\alpha|\times 2n}(\mathbb{R})$
has entries defined in terms of the indices $\\{i\in[2n]\ |\ \alpha_{i}\neq
0\\}=\\{i_{1}<\cdots<i_{r}\\}$ associated with non-zero entries of $\alpha$,
that is,

$\displaystyle(J_{\alpha})_{j,k}=\begin{cases}\delta_{i_{j},k}\qquad&\textrm{
if }j\leq r\\\ 0&\textrm{ otherwise }\ .\end{cases}$ (56)

In other words, $\left(J_{\alpha}\right)_{j,k}=1$ if and only if $k$ is the
position of the $j$-th nonzero element of $\alpha$.

As argued in [60], expressions (49) and (50) determine the inner product
$\langle\Phi_{1},\Phi_{2}\rangle$ and an expression of the form
$\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle$ entirely in terms of covariance
matrices, assuming that the remaining two overlaps
$\langle\Phi_{0},\Phi_{1}\rangle$, $\langle\Phi_{2},\Phi_{0}\rangle$ with a
Gaussian reference state $\Phi_{0}$ are given and non-zero. In this situation,
these quantities can be computed in time $O(n^{3})$.

###  2.7 Gaussian evolution and occupation number measurement

Underlying the known classical simulation algorithms for fermionic linear
optics is the fact that Gaussian unitaries and occupation number measurements
preserve Gaussianity. Explicitly, this can be described as follows: Given a
Gaussian state $\Psi$ with covariance matrix $\Gamma(\Psi)$

  1. (i)

the covariance matrix $\Gamma(U_{R}\Psi)$ of $\Psi$ evolved under the Gaussian
unitary $U_{R}$, $R\in O(2n)$, is given by
$\Gamma(U_{R}\Psi)=R\Gamma(\Psi)R^{T}$.

  2. (ii)

measurement of the observable $a_{j}^{\dagger}a_{j}$ for $j\in[n]$ gives the
outcome $s\in\\{0,1\\}$ with probability

$\displaystyle P_{j}(s)$ $\displaystyle=\|\Pi_{j}(s)\Psi\|^{2}$ (57)
$\displaystyle=\langle\Psi,\Pi_{j}(s)\Psi\rangle$ (58)
$\displaystyle=\frac{1}{2}\left(1+(-1)^{s}\langle\Psi,ic_{2j-1}c_{2j}\Psi\rangle\right)$
(59) $\displaystyle=\frac{1}{2}(1+(-1)^{s}\Gamma_{2j-1,2j})\ .$ (60)

where the third identity follows from the definition

$\displaystyle\Pi_{j}(s)=\frac{1}{2}(1+(-1)^{s}ic_{2j-1}c_{2j})$ (61)

and the last identity follows from the definition of the covariance matrix in
Eq. (43). Given that the measurement outcome is $s\in\\{0,1\\}$, the post-
measurement state

$\displaystyle\Psi(s)$
$\displaystyle=\left(\Pi_{j}\left(s\right)\ket{\Psi}\right)/\sqrt{P_{j}\left(s\right)}$
(62)

is Gaussian with covariance matrix $\Gamma(\Psi(s))$ defined by the lower-
diagonal entries (see e.g., [3, Proof of Lemma 4])

$\displaystyle\Gamma(\Psi(s))_{k,\ell}=\begin{cases}(-1)^{s}&\text{if
}(k,\ell)=(2j,2j-1)\\\
\Gamma_{k,\ell}+\frac{(-1)^{s}}{2P_{j}\left(s\right)}\left(\Gamma_{2j-1,\ell}\Gamma_{2j,k}-\Gamma_{2j-1,k}\Gamma_{2j,\ell}\right)&\text{otherwise}\end{cases}$
(63)

for $k>\ell$.

In particular, the corresponding resulting covariance matrices can be computed
in time $O(n^{3})$ [8, 9, 10, 54] and $O(n^{2})$ [3] for unitary evolution and
measurement, respectively.

###  2.8 The tensor product of two fermionic states

Two density operators $\rho_{j}\in\mathcal{B}(\mathcal{H}^{n_{j}})$,
$j\in[2]$, have a joint extension if and only if there is an element
$\rho\in\mathcal{B}(\mathcal{H}^{n_{1}+n_{2}})$ such that

$\displaystyle\mathop{\mathrm{tr}}\nolimits(c(\alpha_{1}\|\alpha_{2})\rho)$
$\displaystyle=\mathop{\mathrm{tr}}\nolimits(c(\alpha_{1})\rho_{1})\mathop{\mathrm{tr}}\nolimits(c(\alpha_{2})\rho_{2})\qquad\textrm{
for all }\qquad\alpha_{j}\in\\{0,1\\}^{2n_{j}},j\in[2]\ .$ (64)

Here $\alpha_{1}\|\alpha_{2}\in\\{0,1\\}^{2(n_{1}+n_{2})}$ denotes the
concatenation of $\alpha_{1}$ and $\alpha_{2}$. Theorem 1 in [61] implies that
if either $\rho_{1}$ or $\rho_{2}$ is even, then a unique joint extension
$\rho$ of $(\rho_{1},\rho_{2})$ exists. Furthermore, this extension is even if
and only if both $\rho_{1}$ and $\rho_{2}$ are even. Theorem 2 in [61] shows
that if $\rho$ is even and $\rho_{1}$ and $\rho_{2}$ are pure, then $\rho$ is
also pure.

In particular, this means that for states $\Psi_{1},\Psi_{2}$ of definite
parity, there is a unique joint pure extension $\rho=|\Psi\rangle\langle\Psi|$
of
$\left(|\Psi_{1}\rangle\langle\Psi_{1}|,|\Psi_{2}\rangle\langle\Psi_{2}|\right)$.
Since $\rho$ is pure, this also means that $\Psi$ is of definite parity. We
will write $\Psi=\Psi_{1}\tilde{\otimes}\Psi_{2}$ for this state in the
following, and we call $\tilde{\otimes}$ the fermionic tensor product. Note
that $\Psi$ is only defined up to a global phase. It follows immediately from
these definitions that

$\displaystyle\left|\langle x,y|\Psi_{1}\tilde{\otimes}\Psi_{2}\rangle\right|$
$\displaystyle=\left|\langle x,\Psi_{1}\rangle\cdot\langle
y,\Psi_{2}\rangle\right|\qquad\textrm{ for all }\qquad
x\in\\{0,1\\}^{n_{1}}\textrm{ and }y\in\\{0,1\\}^{n_{2}}\ .$ (65)

###### Proof.

Let $x\in\\{0,1\\}^{n_{1}}$ and $y\in\\{0,1\\}^{n_{2}}$ be arbitrary. By
definition, we have

$\displaystyle|x,y\rangle\langle x,y|$
$\displaystyle=\left(\prod_{j=1}^{n_{1}}\frac{1}{2}(I+(-1)^{x_{j}}ic_{2j-1}c_{2j})\right)\left(\prod_{k=1}^{n_{2}}\frac{1}{2}(I+(-1)^{y_{k}}ic_{2n_{1}+2k-1}c_{2n_{1}+2k})\right)$
(66)
$\displaystyle=\left(\sum_{\alpha\in\\{0,1\\}^{2n_{1}}_{+}}\gamma_{x}(\alpha)c(\alpha||0^{2n_{2}})\right)\left(\sum_{\beta\in\\{0,1\\}^{2n_{2}}_{+}}\gamma_{y}(\beta)c(0^{2n_{1}}||\beta)\right)\
.$ (67)

for certain coefficients
$\\{\gamma_{x}(\alpha)\\}_{\alpha\in\\{0,1\\}^{2n_{1}}_{+}}$ and
$\\{\gamma_{y}(\beta)\\}_{\beta\in\\{0,1\\}^{2n_{2}}_{+}}$. Since
$\rho=|\Psi_{1}\tilde{\otimes}\Psi_{2}\rangle\langle\Psi_{1}\tilde{\otimes}\Psi_{2}|$
is an extension of
$(|\Psi_{1}\rangle\langle\Psi_{1}|,|\Psi_{2}\rangle\langle\Psi_{2}|)$ and
$c(\alpha||0^{2n_{2}})c(0^{2n_{1}}||\beta)=i^{-|\alpha|\cdot|\beta|}c(\alpha||\beta)=c(\alpha||\beta)$
for (even-weight) $\alpha\in\\{0,1\\}^{2n_{1}}_{+}$ and
$\beta\in\\{0,1\\}^{2n_{2}}_{+}$, it follows that

$\displaystyle\left|\langle
x,y|\Psi_{1}\tilde{\otimes}\Psi_{2}\rangle\right|^{2}$
$\displaystyle=\sum_{\alpha\in\\{0,1\\}^{2n_{1}}_{+}}\gamma_{x}(\alpha)\sum_{\beta\in\\{0,1\\}^{2n_{2}}_{+}}\gamma_{y}(\beta)\mathop{\mathrm{tr}}\nolimits(c(\alpha\|\beta)\rho)$
(68)
$\displaystyle=\langle\Psi_{1},\left(\sum_{\alpha\in\\{0,1\\}^{2n_{1}}_{+}}\gamma_{x}(\alpha)c(\alpha)\right)\Psi_{1}\rangle\cdot\langle\Psi_{2},\left(\sum_{\beta\in\\{0,1\\}^{2n_{2}}_{+}}\gamma_{y}(\beta)c(\beta)\right)\Psi_{2}\rangle$
(69) $\displaystyle=|\langle x,\Psi_{1}\rangle|^{2}\cdot|\langle
y,\Psi_{2}\rangle|^{2}\ .$ (70)

∎

Refining Eq. (65), (relative) phase information between these matrix elements
can be obtained from the explicit construction of
$\Psi_{1}\tilde{\otimes}\Psi_{2}$ given in [61, Section 3.1] (see also [51,
Proof of Theorem 1]): Consider the isometry

$\displaystyle\begin{matrix}U:&\mathcal{H}^{n_{1}+n_{2}}&\rightarrow&\mathcal{H}^{n_{1}}\otimes\mathcal{H}^{n_{2}}\\\
&\ket{x_{1},\ldots,x_{n_{1}+n_{2}}}&\mapsto&\ket{x_{1},\ldots,x_{m}}\otimes\ket{x_{n_{1}+1},\ldots,x_{n_{1}+n_{2}}}\end{matrix}$
(71)

whose action is given by

$\displaystyle Ua_{j}U^{\dagger}$ $\displaystyle=\begin{cases}a_{j}\otimes
I\qquad&\textrm{ if }j\in[n_{1}]\\\ P_{1}\otimes a_{j-n_{1}}\qquad&\textrm{ if
}j\in\\{n_{1}+1,\ldots,n_{1}+n_{2}\\}\ .\end{cases}$ (72)

where $P_{1}$, the parity operator acting on $\mathcal{H}_{n_{1}}$, introduces
phases. Then

$\displaystyle\Psi_{1}\tilde{\otimes}\Psi_{2}$
$\displaystyle=U^{\dagger}(\Psi_{1}\otimes\Psi_{2})\ .$ (73)

It is straightforward from this definition to check that
$\Psi_{1}\tilde{\otimes}\Psi_{2}$ is the extension of $(\Psi_{1},\Psi_{2})$
and

$\displaystyle\langle x,y|\Psi_{1}\tilde{\otimes}\Psi_{2}\rangle$
$\displaystyle=(-1)^{|x|}\cdot\langle x,\Psi_{1}\rangle\cdot\langle
y,\Psi_{2}\rangle\qquad\textrm{ for all }\qquad x\in\\{0,1\\}^{n_{1}}\textrm{
and }y\in\\{0,1\\}^{n_{2}}\ .$ (74)

We note that the fermionic tensor product preserves Gaussianity in the
following sense.

######  Lemma 2.4.

Let $\Psi_{j}\in\mathcal{G}_{n_{j}}^{+}$ be positive-parity fermionic Gaussian
states for $j\in[2]$. Then
$\Psi_{1}\tilde{\otimes}\Psi_{2}\in\mathcal{G}_{n_{1}+n_{2}}^{+}$, i.e., it is
an even fermionic Gaussian state.

###### Proof.

By definition of an extension (see Eq. (64)) and Wick’s theorem (Eq. (44)),
the tensor product $\Psi=\Psi_{1}\tilde{\otimes}\Psi_{2}$ satisfies

$\displaystyle\langle\Psi,c(\alpha_{1}\|\alpha_{2})\Psi\rangle$
$\displaystyle=\begin{cases}\mathsf{Pf}(\Gamma_{1}[\alpha_{1}])\mathsf{Pf}(\Gamma_{1}[\alpha_{1}])\qquad&\textrm{
if both }|\alpha_{1}|\textrm{ and }|\alpha_{2}|\textrm{ are even }\\\
0&\textrm{ otherwise }\end{cases}$ (75)

for all $\alpha_{j}\in\\{0,1\\}^{2n_{j}}$, where $\Gamma_{j}$ is the
covariance matrix of $\Psi_{j}$ for $j\in[2]$. Because the Pfaffian satisfies

$\displaystyle\mathsf{Pf}(A_{1}\oplus A_{2})$
$\displaystyle=\mathsf{Pf}\begin{pmatrix}A_{1}&0\\\
0&A_{2}\end{pmatrix}=\mathsf{Pf}(A_{1})\mathsf{Pf}(A_{2})$ (76)

for block-matrices, it follows from (75) that the tensor product $\Psi$
satisfies Wick’s theorem (44) with covariance matrix
$\Gamma_{1}\oplus\Gamma_{2}$. In particular, it is Gaussian. ∎

##  3 Tracking relative phases in fermionic linear optics

The covariance matrix $\Gamma(\Psi)$ of a fermionic Gaussian state
$\ket{\Psi}=e^{i\theta}U_{R}\ket{0_{F}}\in\mathcal{G}_{n}$ fully determines
expectation values by Wick’s theorem, which is why Gaussian states and
dynamics are amenable to efficient classical simulation (see Section 2.7).
However, the description of $\Psi$ in terms of the covariance matrix
$\Gamma(\Psi)$ does not capture information on the global phase $e^{i\theta}$
of the state. For processes involving non-Gaussian states expressed as
superpositions of Gaussian states, such phase information needs to be
available for computing norms, expectation values and overlaps.

Here we provide an extended (classical) description of fermionic Gaussian
states that incorporates phase information. A central feature of our
construction is the fact that this description permits to compute overlaps
(including relative phases, i.e., not only the absolute value) of Gaussian
states in an efficient manner.

Our construction is motivated by and relies on expression (49), which relates
the inner product $\langle\Psi_{1},\Psi_{2}\rangle$ of two Gaussian states
$\Psi_{1},\Psi_{2}\in\mathcal{G}_{n}$ to their inner products
$\langle\Psi_{0},\Psi_{1}\rangle$, $\langle\Psi_{0},\Psi_{2}\rangle$ with a
Gaussian reference state $\Psi_{0}\in\mathcal{G}_{n}$ and their covariance
matrices $\Gamma_{0},\Gamma_{1},\Gamma_{2}$. This suggests fixing a reference
state $\Psi_{0}\in\mathcal{G}_{n}$ and using the pair
$(\Gamma(\Psi),\langle\Psi_{0},\Psi\rangle)$ as a classical description of any
state $\ket{\Psi}\in\mathcal{G}_{n}$ relevant in the computation. The problem
with this idea is that $\langle\Psi_{0},\Psi\rangle$ may vanish, preventing
the application of (49). To avoid this problem, we use – instead of a single
state $\Psi_{0}$ – a (potentially) different reference state for each state
$\Psi$. Specifically, we will show that using number states, i.e., states of
the form (18), suffices. This motivates the following definition.

######  Definition 3.1.

Let $\ket{\Psi}=e^{i\theta}U_{R}\ket{0_{F}}\in\mathcal{G}_{n}$ be a Gaussian
state. We call a tuple

$\displaystyle d=(\Gamma,x,r)\in\mathsf{Mat}_{2n\times
2n}(\mathbb{R})\times\\{0,1\\}^{n}\times\mathbb{C}$ (77)

a (valid) description of $\ket{\Psi}$ if the following three conditions hold:

  1. (i)

$\Gamma=\Gamma(\Psi)$ is the covariance matrix of $\ket{\Psi}$.

  2. (ii)

$x\in\\{0,1\\}^{n}$ is such that $\langle x,\Psi\rangle\neq 0$, where
$\ket{x}$ is the number state defined by Eq. (18).

In our algorithms we will in fact ensure that $|\langle x,\Psi\rangle|^{2}\geq
2^{-n}$, i.e., only a subset of valid descriptions is used. A description
$d=(\Gamma,x,r)$ with this property, i.e., satisfying $|r|^{2}\geq 2^{-n}$,
will be called a good description. The restriction to good descriptions is
necessary to make our algorithms work with finite-precision arithmetic.

  3. (iii)

$r=\langle x,\Psi\rangle$.

More explicitly, necessary and sufficient conditions for
$d=(\Gamma(\Psi),x,r)$ to constitute a description of $\Psi$ are that

$\displaystyle r\neq 0\qquad\textrm{ and
}\qquad|r|^{4}=2^{-2n}\mathsf{Det}(\Gamma(\ket{x})+\Gamma(\Psi))$ (78)

because of formula (45) for the overlap of two states and because
$\mathsf{Det}(\cdot)=\mathsf{Pf}^{2}(\cdot)$. Here

$\displaystyle\Gamma(\ket{x})=\bigoplus_{j=1}^{n}\begin{pmatrix}0&(-1)^{x_{j}}\\\
-(-1)^{x_{j}}&0\end{pmatrix}$ (79)

is the covariance matrix of $\ket{x}$. Since a Gaussian state $\Psi$ generally
has non-zero overlap with more than a single occupation number state
$\ket{x}$, there are several distinct valid descriptions of $\Psi$. We will
denote the set of descriptions of $\ket{\Psi}\in\mathcal{G}_{n}$ by
$\mathsf{Desc}(\Psi)$.

We note that a description $d=(\Gamma,x,r)$ uniquely fixes a Gaussian state
$\Psi\in\mathcal{G}_{n}$: The covariance matrix $\Gamma$ determines all
expectation values, and the global phase of $\Psi$ is fixed by the overlap
$\langle x,\Psi\rangle$, i.e., by $r$. Denoting by
$\mathsf{Desc}_{n}=\bigcup_{\Psi\in\mathcal{G}_{n}}\mathsf{Desc}(\Psi)$ the
set of all descriptions of fermionic Gaussian $n$-mode states, this means that
we have a function

$\displaystyle\begin{matrix}\Psi&:\mathsf{Desc}_{n}&\rightarrow&\mathcal{G}_{n}\\\
&d&\mapsto&\Psi(d)\ .\end{matrix}$ (80)

The main result of this section shows that expectation values, overlaps, and
matrix elements of (Majorana) operators with respect to Gaussian states can be
efficiently computed from their classical descriptions. Furthermore, when
evolving a Gaussian state under a Gaussian unitary, the description of the
resulting state can be computed efficiently. The same is true for the post-
measurement state when applying an occupation number measurement.

For evolution, we note that it suffices to consider Gaussian unitaries of the
form $U_{R}$ where $R\in O(2n)$ belongs to the set of generators
$\mathsf{Gen}(O(2n))$ introduced in Section 2.2, that is,

$\displaystyle\mathsf{Gen}(O(2n))=\left\\{R_{j,k}(\vartheta)\ |\
j<k\in[2n],\vartheta\in[0,2\pi)\right\\}\cup\\{R_{j}\\}_{j=1}^{2n}\ .$ (81)

Here $R_{j,k}(\vartheta)$ is a Givens rotation and
$R_{j}=-\mathsf{diag}(\\{(-1)^{\delta_{j,k}}\\}_{k=1}^{2n})$ a reflection. We
note that each element of $\mathsf{Gen}(O(2n))$ can be specified by a tuple
$(j,k,\vartheta)\in[2n]\times[2n]\times[0,2\pi)$ or an index $j\in[2n]$,
respectively. We assume that this parameterization is used in the following
algorithms (but leave this implicit).

To state the properties of our (deterministic) algorithms, it is convenient to
express these as functions.

######  Theorem 3.2 (Overlap, evolution, and measurement).

Let $\Psi(d)\in\mathcal{G}_{n}$ be the Gaussian state associated with a
description $d\in\mathsf{Desc}_{n}$, see Eq. (80). Then the following holds:

  1. (i)

The algorithm
$\mathsf{overlap}:\mathsf{Desc}_{n}\times\mathsf{Desc}_{n}\rightarrow\mathbb{C}$
given in Fig. 7 has runtime $O(n^{3})$ and satisfies

$\displaystyle\mathsf{overlap}(d_{1},d_{2})$
$\displaystyle=\langle\Psi(d_{1}),\Psi(d_{2})\rangle\qquad\textrm{ for all
}\qquad d_{1},d_{2}\in\mathsf{Desc}_{n}\ .$ (82)

  2. (ii)

The algorithm
$\mathsf{evolve}:\mathsf{Desc}_{n}\times\mathsf{Gen}(O(2n))\rightarrow\mathsf{Desc}_{n}$
given in Fig. 9 has runtime $O(n^{3})$ and satisfies

$\displaystyle\Psi(\mathsf{evolve}(d,R))=U_{R}\Psi(d)\qquad\textrm{ for all
}\qquad d\in\mathsf{Desc}_{n}\textrm{ and }R\in\mathsf{Gen}(O(2n))\ ,$ (83)

where $U_{R}$ denotes the Gaussian unitary associated with $R\in O(2n)$.

  3. (iii)

The algorithm
$\mathsf{measureprob}:\mathsf{Desc}_{n}\times[n]\times\\{0,1\\}\rightarrow\mathbb{R}$
given in Fig. 11 has runtime $O(1)$ and satisfies

$\displaystyle\mathsf{measureprob}(d,j,s)=\|\Pi_{j}(s)\Psi(d)\|^{2}\qquad\textrm{
for all }\qquad d\in\mathsf{Desc}_{n},j\in[n],s\in\\{0,1\\}\ ,$ (84)

where $\Pi_{j}(s)=\frac{1}{2}(I+(-1)^{s}ic_{2j-1}c_{2j})$ is the projection
onto the eigenspace of $a_{j}^{\dagger}a_{j}$ with eigenvalue $s$.

  4. (iv)

The algorithm
$\mathsf{postmeasure}:\mathsf{Desc}_{n}\times[n]\times\\{0,1\\}\times[0,1]\rightarrow\mathsf{Desc}_{n}$
given in Fig. 12 has runtime $O(n^{3})$. The algorithm satisfies

$\displaystyle\Psi(\mathsf{postmeasure}(d,j,s,p(d,j,s)))$
$\displaystyle=\frac{\Pi_{j}(s)\Psi(d)}{\|\Pi_{j}(s)\Psi(d)\|}$ (85)

for all $d\in\mathsf{Desc}_{n}$, $j\in[n]$, $s\in\\{0,1\\}$, with
$p(d,j,s)=\|\Pi_{j}(s)\Psi(d)\|^{2}$.

The output of both $\mathsf{evolve}$ and $\mathsf{postmeasure}$ is a good
description for any input.

We argue that descriptions of relevant initial states can be obtained
efficiently. Clearly, this is the case for any state of the form
$\ket{\Psi}=U_{R_{L}}\cdots U_{R_{1}}\ket{0_{F}}$ obtained by applying a
sequence $\\{R_{j}\\}_{j\in[L]}\subset\mathsf{Gen}(O(2n))$ of generators to
the vacuum state $\ket{0_{F}}$: Here we can use the algorithm
$\mathsf{evolve}$ $L$ times, producing a description of $\ket{\Psi}$ in time
$O(Ln^{3})$.

We we will at times need a description of a state $\ket{\Psi}$ but do not
require fixing its global phase. This is the case for example when subsequent
computational states only involve phase-insensitive expressions, e.g., terms
of the form $|\langle\Psi,\Phi\rangle|^{2}$. Such a description can be found
efficiently from the covariance matrix $\Gamma$ of $\ket{\Psi}$. Since the
phase can be fixed arbitrarily, the problem here is to find
$x\in\\{0,1\\}^{n}$ such that $\langle x,\Psi\rangle$ is non-zero.

######  Theorem 3.3.

There is an algorithm $\mathsf{describe}:\mathsf{Mat}_{2n\times
2n}(\mathbb{R})\rightarrow\mathsf{Desc}_{n}$ with runtime $O(n^{3})$ such that
for any covariance matrix $\Gamma$, the state
$\Psi(\mathsf{describe}(\Gamma))$ is a Gaussian state with covariance matrix
$\Gamma$, and $\mathsf{describe}(\Gamma)$ is a good description.

For example, consider states of the form
$\ket{\Phi(\pi,y)}=U_{R_{\pi}}\ket{y}$, where $R_{\pi}\in O(2n)$ is a
permutation matrix specified by an element $\pi\in S_{2n}$ and
$y\in\\{0,1\\}^{n}$. (Such states are used in Ref. [60] to give a fast norm
estimation algorithm, see Section 4.3.) The covariance matrix of this state is
$\Gamma(\pi,y)=R_{\pi}\Gamma(\ket{y})R_{\pi}^{T}$ (with $\Gamma(\ket{y})$
defined by Eq. (79)). We thus conclude that
$\ket{\Psi(\mathsf{describe}(\Gamma(\pi,y)))}$ is proportional to
$\ket{\Phi(\pi_{,}y)}$ with a global phase $e^{i\theta}$ possibly depending on
the pair $(\pi,y)$.

The remainder of this section is devoted to the proofs of Theorem 3.2 and
Theorem 3.3: We describe the algorithms
$\mathsf{evolve},\mathsf{overlap},\mathsf{measureprob},\mathsf{postmeasure}$
and $\mathsf{describe}$ in detail, providing pseudocode, and verify that these
satisfy the desired properties.

###  3.1 Subroutines

Our algorithms make use of subroutines called $\mathsf{findsupport}$,
$\mathsf{relatebasiselements}$, $\mathsf{overlaptriple}$ and
$\mathsf{convert}$ which we describe here.

1:$\Gamma\in\mathsf{Mat}_{2n\times 2n}(\mathbb{R})$ covariance matrix of a
pure Gaussian state  2:function $\mathsf{findsupport}$($\Gamma$)  3:
$\Gamma^{(0)}\leftarrow\Gamma$ 4: $x\leftarrow 0^{n}\in\\{0,1\\}^{n}$.  5: for
$j\leftarrow 1$ to $n$ do $\triangleright$ simulate a measurement of
$a_{j}^{\dagger}a_{j}$ 6:
$q_{j}\leftarrow\frac{1}{2}(I+\Gamma^{(j-1)}_{2j-1,2j})$ 7: if $q_{j}\geq 1/2$
then $\triangleright$ choose the higher-probability outcome  8:
$x_{j}\leftarrow 0$ 9: $p_{j}\leftarrow q_{j}$ 10: else 11: $x_{j}\leftarrow
1$ 12: $p_{j}\leftarrow 1-q_{j}$ 13: $\Gamma^{(j)}\leftarrow
0\in\mathsf{Mat}_{2n\times 2n}(\mathbb{R})$ $\triangleright$ covariance matrix
of the post-measurement state  14:
$\Gamma^{(j)}_{2j,2j-1}\leftarrow(-1)^{x_{j}}$ 15: for $\ell\leftarrow 1$ to
$n-1$ do 16: for $k\leftarrow\ell+1$ to $n$ do 17: if $(k,\ell)\neq(2j,2j-1)$
then 18:
$\Gamma^{(j)}_{k,\ell}\leftarrow\Gamma^{(j-1)}_{k,\ell}+\frac{(-1)^{x_{j}}}{2p_{j}}(\Gamma^{(j-1)}_{2j-1,\ell}\Gamma^{(j-1)}_{2j,k}-\Gamma^{(j-1)}_{2j-1,k}\Gamma^{(j-1)}_{2j,\ell})$
19: $\Gamma^{(j)}\leftarrow\Gamma^{(j)}-(\Gamma^{(j)})^{T}$ 20: return $x$
Figure 1: The algorithm $\mathsf{findsupport}$: Given the covariance matrix
$\Gamma$ of a Gaussian state $\ket{\Psi}$, it computes $x\in\\{0,1\\}^{n}$
such that $\langle x,\Psi\rangle\neq 0$.

The subroutine $\mathsf{findsupport}$ takes as input the covariance matrix
$\Gamma$ of a Gaussian state $\Psi$ and produces a string $x\in\\{0,1\\}^{2n}$
with the property that $\langle x,\Psi\rangle\neq 0$. It is given in Fig. 1.
It has the following properties:

######  Lemma 3.4.

The algorithm $\mathsf{findsupport}:\mathsf{Mat}_{2n\times
2n}(\mathbb{R})\rightarrow\\{0,1\\}^{n}$ runs in time $O(n^{3})$. It satisfies

$\displaystyle|\langle\mathsf{findsupport}(\Gamma(\Psi)),\Psi\rangle|^{2}\geq
2^{-n}\qquad\textrm{ for every }\qquad\Psi\in\mathcal{G}_{n}\ ,$ (86)

where $\Gamma(\Psi)$ is the covariance matrix of $\Psi$.

###### Proof.

The main idea of the algorithm is to mimic a measurement in the number state
basis executed in a sequential manner. Consider the following process: Suppose
we start with the state $\Psi^{(0)}=\Psi$, and then measure
$a_{j}^{\dagger}a_{j}$ successively for $j=1,\dots,n$. Let
$P(x_{j}|x_{1}\cdots x_{j-1})$ denote the conditional probability of observing
the outcome $x_{j}\in\\{0,1\\}$ (when measuring $a_{j}^{\dagger}a_{j})$, given
that the previous measurements yielded $(x_{1},\ldots,x_{j-1})$. According to
Born’s rule, this is given by

$\displaystyle P(x_{j}|x_{1}\cdots x_{j-1})$
$\displaystyle=\langle\Psi^{(j-1)}_{x_{1}\cdots
x_{j-1}},\Pi_{j}(x_{j})\Psi^{(j-1)}_{x_{1}\cdots x_{j-1}}\rangle$ (87)

where $\Psi^{(j-1)}_{x_{1}\cdots x_{j-1}}$ is the post-measurement state after
the first $(j-1)$ measurements. The probability of observing the sequence
$x\in\\{0,1\\}^{n}$ of outcomes then is

$\displaystyle|\langle x,\Psi\rangle|^{2}$
$\displaystyle=\prod_{j=1}^{n}P(x_{j}|x_{1}\cdots x_{j-1})$ (88)

by Bayes’ rule.

The algorithm $\mathsf{findsupport}$ simulates this process: For each
$j\in[n]$, the quantity $q_{j}$ computed in line 6 is equal to the conditional
probability $P(0|x_{1}\cdots x_{j-1})$ that the $j$-th measurement results in
the outcome $0$. Lines 7–12 ensure that the outcome $x_{j}\in\\{0,1\\}$ with
higher probability of occurrence is selected at each step, guaranteeing
Property (86) (because of Eq. (88)). The matrix $\Gamma^{(j)}$ computed in
steps 13–19 is the covariance matrix of the post-measurement state
$\Psi^{(j)}_{x_{1}\cdots x_{j}}$.

Each measurement is thus realized in time $O(n^{2})$ yielding the overall
complexity of $O(n^{3})$. ∎

The algorithm $\mathsf{relatebasiselements}$ is more straightforward: Given
$x,y\in\\{0,1\\}^{n}$, it outputs
$(\alpha,\vartheta)\in\\{0,1\\}^{2n}\times\mathbb{R}$ such that
$c(\alpha)\ket{x}=e^{i\vartheta}\ket{y}$. That is, it finds a Majorana
monomial $c(\alpha)$ which maps the basis state $\ket{x}$ to $\ket{y}$ up to a
phase and computes the corresponding phase. In Fig. 2 we give pseudocode for
this algorithm.

1:$x\in\\{0,1\\}^{n}$, $y\in\\{0,1\\}^{n}$ 2:function
$\mathsf{relatebasiselements}$($x$,$y$)  3: $\alpha\leftarrow
0^{2n}\in\\{0,1\\}^{2n}$.  4: for $j\leftarrow 1$ to $n$ do 5:
$\alpha_{2j-1}\leftarrow x_{j}\oplus y_{j}$ 6:
$\vartheta\leftarrow\frac{\pi}{4}|x\oplus y|+\pi\sum_{j=1}^{n}(x\oplus
y)_{j}\eta_{j}(x)$ 7: return $(\alpha,\vartheta)$ $\triangleright$
$(\alpha,\vartheta)$ is such that $c(\alpha)\ket{x}=e^{i\vartheta}\ket{y}$
Figure 2: Given $x,y\in\\{0,1\\}^{n}$, the algorithm
$\mathsf{relatebasiselements}$ computes $\alpha\in\\{0,1\\}^{2n}$ and
$\vartheta\in\mathbb{R}$ such that $c(\alpha)\ket{x}=e^{i\vartheta}\ket{y}$.
The expression $\eta_{j}(x)$ in line 6 is defined by Eq. (20), and $x\oplus
y\in\\{0,1\\}^{n}$ denotes the bitwise addition modulo two for
$x,y\in\\{0,1\\}^{n}$.

######  Lemma 3.5.

The algorithm
$\mathsf{relatebasiselements}:\\{0,1\\}^{n}\rightarrow\\{0,1\\}^{2n}\times\mathbb{C}$
runs in time $O(n)$ and satisfies

$\displaystyle c(\alpha)\ket{x}=e^{i\vartheta}\ket{y}\textrm{ where
}(\alpha,\vartheta)=\mathsf{relatebasiselements}(x,y)\,\quad\textrm{ for all
}\quad x,y\in\\{0,1\\}^{n}\ .$ (89)

###### Proof.

Let $x,y\in\\{0,1\\}^{n}$ be arbitrary. Define

$\displaystyle\alpha_{2j-1}$ $\displaystyle=x_{j}\oplus y_{j}\qquad\textrm{
and }\qquad\alpha_{2j}=0\qquad\textrm{ for }\qquad j\in[n]\ ,$ (90)

as in line 5 of algorithm $\mathsf{relatebasiselements}$. Then

$\displaystyle c(\alpha)\ket{y}$
$\displaystyle=i^{|\alpha|\cdot(|\alpha|-1)/2}c_{1}^{x_{1}\oplus
y_{1}}c_{3}^{x_{2}\oplus y_{2}}\cdots c_{2n-1}^{x_{n}\oplus y_{n}}\ket{y}$
(91) $\displaystyle=i^{\left(\sum_{j=1}^{n}x_{j}\oplus
y_{j}\right)\cdot\left(\left(\sum_{j=1}^{n}x_{j}\oplus
y_{j}\right)-1\right)/2}(-1)^{\sum_{j=1}^{n}(x_{j}\oplus
y_{j})\eta_{j}(x)}\ket{y\oplus(x\oplus y)}$ (92) $\displaystyle=i^{|x\oplus
y|\cdot(|x\oplus y|-1)/2}(-1)^{\sum_{j=1}^{n}(x\oplus
y)_{j}\eta_{j}(x)}\ket{x}$ (93)

where in the second identity, we used that

$\displaystyle c_{2j-1}\ket{x}$ $\displaystyle=(-1)^{\eta_{j}(x)}\ket{x\oplus
e_{j}}\qquad\textrm{ for all }\qquad x\in\\{0,1\\}^{n}\textrm{ and }j\in[n]$
(94)

because of Eq. (19). Because $i^{|x\oplus y|\cdot(|x\oplus
y|-1)/2}(-1)^{\sum_{j=1}^{n}(x\oplus y)_{j}\eta_{j}(x)}=e^{i\vartheta}$ for

$\displaystyle\vartheta=\frac{\pi}{4}|x\oplus y|\cdot(|x\oplus
y|-1)+\pi\sum_{j=1}^{n}(x\oplus y)_{j}\eta_{j}(x)\ ,$ (95)

comparison with line 6 of the algorithm $\mathsf{relatebasiselements}$ gives
the claim. ∎

The algorithm $\mathsf{overlaptriple}$ takes covariance matrices
$\Gamma_{0},\Gamma_{1},\Gamma_{2}$ of three Gaussian states
$\Phi_{0},\Phi_{1},\Phi_{2}$ of the same parity $\sigma\in\\{-1,1\\}$ and
$\alpha\in\\{0,1\\}_{+}^{n}$, as well as overlaps
$u=\langle\Phi_{0},\Phi_{1}\rangle$,
$v=\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle$ (which both have to be non-zero),
and computes the overlap $\langle\Phi_{2},\Phi_{0}\rangle$. It is obtained by
direct application of the formula (50). For completeness, we include
pseudocode in Fig. 3. Since this algorithm involves computing Pfaffians of
matrices that have size linear in $n$, its runtime is $O(n^{3})$.

1:$\Gamma_{j}$ covariance matrix of a Gaussian state $\Phi_{j}$ for $j=0,1,2$
2:$\sigma=\mathsf{Pf}(\Gamma_{0})=\mathsf{Pf}(\Gamma_{1})=\mathsf{Pf}(\Gamma_{2})$
3:$\alpha\in\\{0,1\\}^{n}_{+}$ such that
$\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle\neq 0$
4:$\langle\Phi_{0},\Phi_{1}\rangle\neq 0$
5:$u=\langle\Phi_{0},\Phi_{1}\rangle$ and
$v=\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle$ 6:function
$\mathsf{overlaptriple}$($\Gamma_{0},\Gamma_{1},\Gamma_{2},\alpha,u,v$)  7:
$R_{\alpha}\in\mathsf{Mat}_{(6n+|\alpha|)\times(6n+|\alpha|)}(\mathbb{R})$ 8:
$R_{\alpha}\leftarrow\textrm{ evaluate Eq.~{}\eqref{eq:Ralphadefinitionm}}$ 9:
$o\leftarrow u^{-1}\cdot v^{-1}\cdot\sigma\cdot
4^{-n}i^{n+|\alpha|\cdot(|\alpha|-1)/2}\mathsf{Pf}(R_{\alpha})$
$\triangleright$ compute the overlap $\langle\Phi_{2},\Phi_{0}\rangle$ 10:
return $o$ Figure 3: The algorithm $\mathsf{overlaptriple}$ takes as input the
covariance matrices $\Gamma_{j}$ of three Gaussian states $\Phi_{j}$,
$j=0,1,2$ with identical parity, $\alpha\in\\{0,1\\}_{+}^{n}$ and the overlaps
$\langle\Phi_{0},\Phi_{1}\rangle$, $\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle$.
The latter have to be non-zero. The algorithm computes the overlap
$\langle\Phi_{2},\Phi_{0}\rangle$ using Eq. (50).

Fig. 4 gives a graphical representation of what the algorithm
$\mathsf{overlaptriple}$ achieves. These graphical representations will be
helpful to construct and analyze other algorithmic building blocks.

(a) The input to the algorithm $\mathsf{overlaptriple}$ consists of
(descriptions of) three Gaussian states $\Phi_{0},\Phi_{1},\Phi_{2}$ and
$\alpha\in\\{0,1\\}^{2n}$, together with overlaps
$u=\langle\Phi_{0},\Phi_{1}\rangle$ and
$v=\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle$ that are both non-zero.

(b) Applying $\mathsf{overlaptriple}$ provides the inner product
$\langle\Phi_{0},\Phi_{2}\rangle$. In this diagrammatic representation, this
completes the triangle with vertices $\Phi_{0},\Phi_{1},\Phi_{2}$.

Figure 4: A graphical representation of the functionality provided by the
algorithm $\mathsf{overlaptriple}$. Solid lines represent inner products that
are given / have been computed, and are non-zero. Inner products of the form
$\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle$ are represented by arrows.

The algorithm $\mathsf{convert}$ takes a description $d=(\Gamma,x,r)$ of a
Gaussian state $\Psi(d)$ and $y\in\\{0,1\\}^{n}$ such that $\langle
y,\Psi(d)\rangle\neq 0$, and outputs a description $d^{\prime}=(\Gamma,y,s)$
of the same state. In other words, it converts a description $d$ of the state
involving the reference state $\ket{x}$ to a description $d^{\prime}$ of the
same state but involving a different reference state $\ket{y}$. In Fig. 5 we
give pseudocode for this algorithm.

1:$d=(\Gamma,x,r)\in\mathsf{Desc}_{n}$, $y\in\\{0,1\\}^{n}$ such that $\langle
y,\Psi(d)\rangle\neq 0$ 2:function $\mathsf{convert}$($d$,$y$)  3:
$(\alpha,\vartheta)\leftarrow\mathsf{relatebasiselements}(y,x)$
$\triangleright$ find $(\alpha,\vartheta)$ such that
$c(\alpha)\ket{y}=e^{i\vartheta}\ket{x}$ 4: $\Gamma_{0}\leftarrow\Gamma$,
$\Gamma_{1}\leftarrow\Gamma(\ket{x})$, $\Gamma_{2}\leftarrow\Gamma(\ket{y})$
$\triangleright$ covariance matrices of $\Psi(d)$, $\ket{x}$ and $\ket{y}$ 5:
$u\leftarrow\overline{r}$ $\triangleright$ $u=\langle\Psi(d),x\rangle$ 6:
$v\leftarrow e^{i\vartheta}$ $\triangleright$ $v=\langle x,c(\alpha)y\rangle$
7:
$w\leftarrow\mathsf{overlaptriple}(\Gamma_{0},\Gamma_{1},\Gamma_{2},\alpha,u,v)$
$\triangleright$ compute the overlap $\langle y,\Psi(d)\rangle$ 8: return
$(\Gamma,y,w)$ Figure 5: The algorithm $\mathsf{convert}$ takes a description
$d\in\mathsf{Desc}_{n}$ and $y\in\\{0,1\\}^{n}$ such that $\langle
y,\Psi(d)\rangle\neq 0$. It outputs a description
$d^{\prime}\in\mathsf{Desc}_{n}$ of $\Psi(d)$ such that the second entry of
$d^{\prime}$ is equal to $y$, i.e., $d^{\prime}=(\Gamma,y,s)$. It makes use of
the subroutines $\mathsf{relatebasiselements}$ and $\mathsf{overlaptriple}$.
For $x\in\\{0,1\\}^{n}$, $\Gamma(\ket{x})$ denotes the covariance matrix of
the state $\ket{x}$, see Eq. (79).

The algorithm $\mathsf{convert}$ is illustrated in Fig. 6.

(a) The input to the algorithm $\mathsf{convert}$ specifies a Gaussian state
$\Psi$, $x\in\\{0,1\\}^{n}$ such that $\langle x,\Psi\rangle\neq 0$, the value
$r=\langle x,\Psi\rangle$ and an element $y\in\\{0,1\\}^{n}$ such that
$\langle y,\Psi\rangle\neq 0$. The value of $\langle y,\Psi\rangle$ is not
given.

(b) The algorithm applies the subroutine $\mathsf{relatebasiselements}$ to
find $(\alpha,\vartheta)$ such that $c(\alpha)\ket{y}=e^{i\vartheta}\ket{x}$.
In particular, after this step, the value $\langle
x,c(\alpha)y\rangle=e^{i\vartheta}$ is known and it is non-zero.

(c) The algorithm then applies the subroutine $\mathsf{overlaptriple}$ to
compute $w=\langle y,\Psi\rangle$. The triple $(\Gamma,y,w)$ is a valid
description of $\Psi$.

Figure 6: An illustration of the algorithm $\mathsf{convert}$. Dotted lines
represent inner products that are non-zero, but that are not provided / have
not yet been computed by the algorithm.

######  Lemma 3.6.

The algorithm
$\mathsf{convert}:\mathsf{Desc}_{n}\times\\{0,1\\}^{n}\rightarrow\mathsf{Desc}_{n}$
given in Fig. 5 runs in time $O(n^{3})$. Assume that $d\in\mathsf{Desc}_{n}$
and $y\in\\{0,1\\}^{n}$ satisfy $\langle y,\Psi(d)\rangle\neq 0$. Then

$\displaystyle\Psi(\mathsf{convert}(d,y))=\Psi(d)\ .$ (96)

Furthermore, denoting the output of $\mathsf{convert}(d,y)$ by
$d^{\prime}=(\Gamma^{\prime},y^{\prime},s^{\prime})$ we have

$\displaystyle y^{\prime}=y$ (97)

as well as

$\displaystyle|s^{\prime}|^{2}=|\langle
y^{\prime},\Psi(d)\rangle|^{2}=|\langle y,\Psi(d)\rangle|^{2}\ .$ (98)

###### Proof.

Let us denote the input to $\mathsf{convert}$ by $(d,y)$, where
$d=(\Gamma,x,r)\in\mathsf{Desc}_{n}$ and $y\in\\{0,1\\}^{n}$. Then

$\displaystyle\langle x,\Psi(d)\rangle\neq 0$ (99)

since $d$ is a description of $\Psi(d)$. Furthermore, for $(\alpha,\vartheta)$
as defined in line 3 we have

$\displaystyle\langle x,c(\alpha)y\rangle=e^{i\vartheta}\neq 0$ (100)

by definition of the algorithm $\mathsf{relatebasiselements}$. In line 4 of
$\mathsf{convert}$, the matrices $\Gamma_{j}$, $j\in[3]$ are the covariance
matrices of the states

$\displaystyle(\Phi_{0},\Phi_{1},\Phi_{2})$
$\displaystyle=(\Psi(d),\ket{x},\ket{y})\ .$ (101)

We note that Eq. (99) and the assumption $\langle y,\Psi(d)\rangle\neq 0$
imply that these three states have identical parity. The value $w$ computed in
line 7 using $\mathsf{overlaptriple}$ is equal to the overlap

$\displaystyle w$ $\displaystyle=\langle\Phi_{2},\Phi_{0}\rangle=\langle
y,\Psi(d)\rangle\ ,$ (102)

because

$\displaystyle u$ $\displaystyle=\overline{r}=\overline{\langle
x,\Psi(d)\rangle}=\langle\Psi(d),x\rangle=\langle\Phi_{0},\Phi_{1}\rangle\neq
0$ (103) $\displaystyle v$ $\displaystyle=e^{i\vartheta}=\langle
x,c(\alpha)y\rangle=\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle\neq 0$ (104)

by Eqs. (99) and (100). Eq. (102) together with the assumption $\langle
y,\Psi(d)\rangle\ \neq 0$ show that the output $(\Gamma,y,w)$ is a description
of $\Psi(d)$. This completes the proof of Eq. (96).

Eq. (98) is trivially satisfied because

$\displaystyle s^{\prime}=w=\langle y,\Psi(d)\rangle\ .$ (105)

The complexity of the algorithm is dominated by $\mathsf{overlaptriple}$,
which takes time $O(n^{3})$. ∎

###  3.2 Computing overlaps and descriptions of evolved/measured states

Based on the subroutines $\mathsf{findsupport}$,
$\mathsf{relatebasiselements}$, $\mathsf{overlaptriple}$ and
$\mathsf{convert}$, we can now describe our main algorithms
$\mathsf{overlap}$, $\mathsf{evolve}$, $\mathsf{measureprob}$ and
$\mathsf{postmeasure}$ for overlaps, Gaussian unitary evolution, to compute
the outcome probability and the post-measurement state when measuring the
occupation number, respectively. We give pseudocode for each algorithm and
establish the associated claims.

We give pseudocode for the algorithm $\mathsf{overlap}$ in Fig. 7 and we
illustrate it in Fig. 8.

1:$d_{1}=(\Gamma_{1},x_{1},r_{1}),d_{2}=(\Gamma_{2},x_{2},r_{2})\in\mathsf{Desc}_{n}$
2:function $\mathsf{overlap}$($d_{1}$,$d_{2}$)  3:
$\sigma_{1}\leftarrow\mathsf{Pf}(\Gamma_{1})$ $\triangleright$ compute the
parity $\sigma_{j}$ of $\Psi(d_{j})$ 4:
$\sigma_{2}\leftarrow\mathsf{Pf}(\Gamma_{2})$ 5: if $\sigma_{1}\neq\sigma_{2}$
then 6: return 0 $\triangleright$ states with different parities have zero
overlap  7:
$(\alpha,\vartheta)\leftarrow\mathsf{relatebasiselements}(x_{2},x_{1})$
$\triangleright$ $(\alpha,\vartheta)$ satisfies
$c(\alpha)\ket{x_{2}}=e^{i\vartheta}\ket{x_{1}}$ 8:
$\Gamma^{\prime}_{0}\leftarrow\Gamma_{1}$,
$\Gamma^{\prime}_{1}\leftarrow\Gamma(\ket{x_{1}})$,
$\Gamma^{\prime}_{2}\leftarrow\Gamma_{2}$ $\triangleright$ covariance matrices
of $\Psi(d_{1}),\ket{x_{1}}$ and $\Psi(d_{2})$ 9:
$u\leftarrow\overline{r}_{1}$ $\triangleright$
$u=\langle\Psi(d_{1}),x_{1}\rangle$ 10: $v\leftarrow e^{i\vartheta}r_{2}$
$\triangleright$ $v=\langle x_{1},c(\alpha)\Psi(d_{2})\rangle$ 11:
$w\leftarrow\mathsf{overlaptriple}(\Gamma^{\prime}_{0},\Gamma^{\prime}_{1},\Gamma^{\prime}_{2},\alpha,u,v)$
$\triangleright$ compute the overlap $\langle\Psi(d_{2}),\Psi(d_{1})\rangle$
12: return $\overline{w}$ Figure 7: The algorithm $\mathsf{overlap}$ takes
descriptions $d_{1},d_{2}\in\mathsf{Desc}_{n}$ and outputs the overlap
$\langle\Psi(d_{1}),\Psi(d_{2})\rangle$.

(a) The input to the algorithm $\mathsf{overlap}$ consists of (descriptions
of) two Gaussian $\Psi_{1},\Psi_{2}$, $x_{1},x_{2}\in\\{0,1\\}^{n}$ and the
overlaps $r_{j}=\langle x_{j},\Psi_{j}\rangle$, $j\in[2]$. The latter are both
assumed to be non-zero.

(b) The algorithm uses the subroutine $\mathsf{relatebasiselements}$ to find
$(\alpha,\vartheta)$ such that
$c(\alpha)\ket{x_{2}}=e^{i\vartheta}\ket{x_{1}}$. In particular, this means
that the value $\langle x_{1},c(\alpha)x_{2}\rangle=e^{i\vartheta}$ is
computed, and it is non-zero. Furthermore, this implies that $\langle
x_{1},c(\alpha)\Psi_{2}\rangle=\langle
c(\alpha)x_{1},\Psi_{2}\rangle=e^{i\vartheta}\langle
x_{2},\Psi_{2}\rangle=e^{i\vartheta}r_{2}$ is also known and non-zero.

(c) In the last step of the algorithm, the subroutine $\mathsf{overlaptriple}$
is applied to complete a triangle: This amounts to computing
$w=\langle\Psi_{2},\Psi_{1}\rangle$. The algorithm returns the complex
conjugate $\overline{w}=\langle\Psi_{1},\Psi_{2}\rangle$.

Figure 8: An illustration of the algorithm $\mathsf{overlap}$.

######  Lemma 3.7.

The algorithm
$\mathsf{overlap}:\mathsf{Desc}_{n}\times\mathsf{Desc}_{n}\rightarrow\mathbb{C}$
given in Fig. 7 runs in time $O(n^{3})$. It satisfies

$\displaystyle\mathsf{overlap}(d_{1},d_{2})$
$\displaystyle=\langle\Psi(d_{1}),\Psi(d_{2})\rangle\qquad\textrm{ for all
}\qquad d_{1},d_{2}\in\mathsf{Desc}_{n}\ .$ (106)

###### Proof.

Let $d_{j}=(\Gamma_{j},x_{j},r_{j})\in\mathsf{Desc}_{n}$ for $j\in[2]$. Then

$\displaystyle r_{j}$ $\displaystyle=\langle x_{j},\Psi(d_{j})\rangle\neq
0\qquad\textrm{ for }\qquad j\in[2]\ ,$ (107)

by assumption.

Line 5 treats the case where $\Psi(d_{1})$ and $\Psi(d_{2})$ have different
parity, and are thus orthogonal. Starting from line 7, we can hence assume
that the parities $\sigma_{1},\sigma_{2}$ of $\Psi(d_{1}),\Psi(d_{2})$ are
identical, $\sigma=\sigma_{1}=\sigma_{2}$. By Eq. (107), this implies that
both $\ket{x_{1}}$ and $\ket{x_{2}}$ also have parity $\sigma$, that is,

$\displaystyle\sigma(\ket{x_{1}})=\sigma(\ket{x_{2}})=\sigma(\ket{\Psi(d_{1})})=\sigma(\ket{\Psi(d_{2})})\
.$ (108)

By definition of $\mathsf{relatebasiselements}$, the pair $(\alpha,\vartheta)$
computed in line 7 satisfies

$\displaystyle c(\alpha)\ket{x_{2}}=e^{i\vartheta}\ket{x_{1}}\ .$ (109)

Consider the triple of states

$\displaystyle(\Phi_{0},\Phi_{1},\Phi_{2})$
$\displaystyle=(\Psi(d_{1}),\ket{x_{1}},\Psi(d_{2}))\ .$ (110)

Then the matrices $\Gamma_{j}^{\prime}$, $j\in[3]$ defined in line 8 of the
algorithm are the covariance matrices of $\Phi_{j}$, $j\in\\{0,1,2\\}$. We
have

$\displaystyle u$ $\displaystyle=\overline{r}_{1}=\overline{\langle
x_{1},\Psi(d_{1})\rangle}=\langle\Psi(d_{1}),x_{1}\rangle=\langle\Phi_{0},\Phi_{1}\rangle\neq
0\ ,$ (111)

by Eq. (107), and similarly

$\displaystyle v$ $\displaystyle=e^{i\vartheta}r_{2}$ (112)
$\displaystyle=e^{i\vartheta}\langle x_{2},\Psi(d_{2})\rangle$ (113)
$\displaystyle=\langle e^{-i\vartheta}x_{2},\Psi(d_{2})\rangle$ (114)
$\displaystyle=\langle c(\alpha)x_{1},\Psi(d_{2})\rangle$ with (109) (115)
$\displaystyle=\langle x_{1},c(\alpha)\Psi(d_{2})\rangle$ because $c(\alpha)$
is self-adjoint (116) $\displaystyle=\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle$
(117) $\displaystyle\neq 0$ because $r_{2}\neq 0$ by Eq. (107) . (118)

Furthermore, by Eq. (108) the states $(\Phi_{0},\Phi_{1},\Phi_{2})$ have
identical parity. It thus follows from the properties of
$\mathsf{overlaptriple}$ that the quantity $w$ computed in step line 11 is
equal to

$\displaystyle w$ $\displaystyle=\langle\Psi(d_{2}),\Psi(d_{1})\rangle\ .$
(119)

Since the output of the algorithm is the complex conjugate
$\overline{w}=\langle\Psi(d_{1}),\Psi(d_{2})\rangle$, this implies the claim
(106).

The runtime of $\mathsf{overlap}$ is dominated by $\mathsf{overlaptriple}$,
hence it is of order $O(n^{3})$. ∎

We give pseudocode for the algorithm $\mathsf{evolve}$ in Fig. 9 and we
illustrate it in Fig. 10(a).

1:$R\in\mathsf{Gen}(O(2n))$ 2:$d=(\Gamma,x,r)\in\mathsf{Desc}_{n}$ 3:function
$\mathsf{evolve}$($R,d$)  4: $\Gamma_{0}\leftarrow R\Gamma R^{T}$
$\triangleright$ covariance matrix of $U_{R}\Psi(d)$ 5:
$y\leftarrow\mathsf{findsupport}(\Gamma_{0})$ $\triangleright$ gives $y$ such
that $|\langle y,U_{R}\Psi(d)\rangle|^{2}\geq 2^{-n}$ 6: if
$R=R_{j,k}(\vartheta)$ for $j<k\in[2n],\vartheta\in\mathbb{R}$ then
$\triangleright$ $R=R_{j,k}(\vartheta)$ is a Givens rotation  7: if
$\cos^{2}(\vartheta/2)\geq 1/2$ then 8: $z\leftarrow x$ 9:
$s\leftarrow\cos(\vartheta/2)$ 10: else 11: $z\leftarrow x\oplus e_{j}\oplus
e_{k}$ 12: $\beta\leftarrow\beta_{j}(x)+\beta_{k}(x)$ 13: $s\leftarrow
e^{i\pi\beta}\sin(\vartheta/2)$ 14: if $R=R_{j}$ for $j\in[2n]$ then
$\triangleright$ $R=R_{j}$ is a reflection  15: $z\leftarrow x\oplus e_{j}$
16: $\beta\leftarrow\beta_{j}(x)$ 17: $s\leftarrow e^{i\pi\beta}$
$\triangleright$ $(s,z)$ satisfies $s=\langle z,U_{R}x\rangle$, $|s|^{2}\geq
1/2$ 18: $(\alpha,\gamma)\leftarrow\mathsf{relatebasiselements}(y,z)$
$\triangleright$ $(\alpha,\gamma)$ satisfies
$c(\alpha)\ket{y}=e^{i\gamma}\ket{z}$ 19: $\Gamma_{1}\leftarrow
R\Gamma(\ket{x})R^{T}$, $\Gamma_{2}\leftarrow\Gamma(\ket{y})$ $\triangleright$
covariance matrices of $U_{R}\ket{x}$ and $\ket{y}$ 20:
$u\leftarrow\overline{r}$ $\triangleright$ $u=\langle U_{R}\Psi,U_{R}x\rangle$
21: $v\leftarrow e^{i\gamma}\overline{s}$ $\triangleright$ $v=\langle
U_{R}x,c(\alpha)y\rangle$ 22:
$w\leftarrow\mathsf{overlaptriple}(\Gamma_{0},\Gamma_{1},\Gamma_{2},\alpha,u,v)$
$\triangleright$ compute the overlap $\langle y,U_{R}\Psi(d)\rangle$ 23:
return $(\Gamma_{0},y,w)$ $\triangleright$ return a description of
$U_{R}\ket{\Psi(d)}$ Figure 9: The algorithm $\mathsf{evolve}$ takes a
description $d\in\mathsf{Desc}_{n}$ and an orthogonal matrix
$R\in\mathsf{Gen}(O(2n))$ associated with the Gaussian unitary $U_{R}$ and
computes a description for the state $U_{R}\Psi(d)$. In this algorithm, the
functions $\beta_{s}:\\{0,1\\}^{n}\rightarrow\mathbb{R}$ for $s\in[n]$ are
defined as
$\beta_{s}(x)=\eta_{s}(x)+\left(x_{s}-\frac{1}{2}\right)\cdot(s+1)$,
$x\in\\{0,1\\}^{n}$ with $\eta_{s}(x)$ given in Eq. (20).

(a) The input to the algorithm $\mathsf{evolve}$ is (the description of) a
Gaussian state $\Psi$, $x\in\\{0,1\\}^{n}$ and the non-zero overlap $r=\langle
x,\Psi\rangle$, as well as an element $R\in\mathsf{Gen}(O(2n))$ associated
with a Gaussian unitary $U_{R}$.

(b) By unitarity of $U_{R}$, the input data also provides the inner product
$\langle U_{R}x,U_{R}\Psi\rangle=\langle x,\Psi\rangle=r$, which is non-zero.

(c) The algorithm invokes the subroutine $\mathsf{findsupport}$ applied to the
covariance matrix $\Gamma_{0}$ of the evolved state $U_{R}\Psi$ in order to
find an element $y\in\\{0,1\\}^{n}$ such that $\langle y,U_{R}\Psi\rangle\neq
0$. (The value of this inner product is not computed/available at this point.)

(d) Using the fact the action of $U_{R}$ for a generator
$R\in\mathsf{Gen}(O(2n))$ is local, the algorithm determines an element
$z\in\\{0,1\\}^{n}$ such that $s=\langle z,U_{R}x\rangle\neq 0$, and computes
the value $s$.

(e) The algorithm then uses the subroutine $\mathsf{relatebasiselements}$ to
find $(\alpha,\gamma)$ such that $c(\alpha)\ket{y}=e^{i\gamma}\ket{z}$. This
means that the inner product $\langle z,c(\alpha)y\rangle=e^{i\gamma}$ is
known and non-zero. Since $\langle z,U_{R}x\rangle$ is known and non-zero (as
ensured by the previous step (10(d))), the value $\langle
U_{R}x,c(\alpha)y\rangle=e^{i\gamma}\overline{\langle z,U_{R}x\rangle}$ is
also non-zero and can be computed.

(f) In the last step, the subroutine $\mathsf{overlaptriple}$ is used to
compute the quantity $w=\langle y,U_{R}\Psi\rangle$. It is non-zero by step
(10(c)). Thus $(\Gamma_{0},y,w)$ is a valid description of $U_{R}\Psi$.

Figure 10: An illustration of the algorithm $\mathsf{evolve}$. Dotted lines
correspond to inner products whose value is non-zero, but has not been
computed at that stage of the algorithm.

######  Lemma 3.8.

The algorithm
$\mathsf{evolve}:\mathsf{Gen}(O(2n))\times\mathsf{Desc}_{n}\rightarrow\mathsf{Desc}_{n}$
given in Fig. 9 runs in time $O(n^{3})$. Consider an arbitrary generator
$R\in\mathsf{Gen}(O(2n))$ and a description $d\in\mathsf{Desc}_{n}$. Then

$\displaystyle\Psi(\mathsf{evolve}(R,d))=U_{R}\Psi(d)\ ,$ (120)

that is, the output of $\mathsf{evolve}$ is a description of the evolved state
$U_{R}\Psi(d)$. Furthermore, denoting the output by
$d^{\prime}=(\Gamma^{\prime},x^{\prime},r^{\prime})=\mathsf{evolve}(R,d)$ we
have

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

###### Proof.

Let us denote the input of $\mathsf{evolve}$ by $(R,d)$ where
$R\in\mathsf{Gen}(O(2n))$ and $d=(\Gamma,x,r)\in\mathsf{Desc}_{n}$. The state
$U_{R}\Psi(d)$ has covariance matrix $\Gamma_{0}=R\Gamma R^{T}$ computed in
line 4 (see Section 2.7). By the properties of $\mathsf{findsupport}$ (see
Lemma 3.4), the state $\ket{y}$ with
$y=\mathsf{findsupport}(\Gamma_{0})\in\\{0,1\\}^{n}$ computed in line 5 is
such that

$\displaystyle|\langle y,U_{R}\Psi(d)\rangle|^{2}\geq 2^{-n}\ .$ (122)

In particular, it is non-zero. The remainder of the algorithm computes
$\langle y,U_{R}\Psi(d)\rangle$.

We first show the following:

######  Claim 3.9.

Lines 6–17 compute $(z,s)\in\\{0,1\\}^{n}\times\mathbb{C}$ such that

$\displaystyle|\langle z,U_{R}x\rangle|^{2}\geq 1/2$ (123)

and

$\displaystyle s=\langle z,U_{R}x\rangle\ .$ (124)

###### Proof.

Here we are using the fact that for any generator $R\in\mathsf{Gen}(O(2n))$,
the associated Gaussian unitary $U_{R}$ has a local action on the mode
operators. In particular, we can easily compute the image $U_{R}\ket{x}$ of a
number state $\ket{x}$ under $U_{R}$. We distinguish two cases:

  1. (i)

$R=R_{j,k}(\vartheta)$, $j<k\in[2n]$, $\vartheta\in[0,2\pi)$ is a Givens-
rotation (see Lines 6–13): In this case, $R$ is associated with the unitary
evolution operator

$\displaystyle
U_{j,k}=\exp(\vartheta/2c_{j}c_{k})=\cos(\vartheta/2)I+\sin(\vartheta/2)c_{j}c_{k}\
.$ (125)

It maps a basis state $\ket{x}$, $x\in\\{0,1\\}^{n}$, to

$\displaystyle
U_{j,k}(\vartheta)\ket{x}=\cos(\vartheta/2)\ket{x}+e^{i\pi(\beta_{j}(x)+\beta_{k}(x))}\sin(\vartheta/2)\ket{x\oplus
e_{j}\oplus e_{k}}\ $ (126)

where we introduced the quantities

$\displaystyle\beta_{s}(x)=\eta_{s}(x)+\left(x_{s}-\frac{1}{2}\right)\cdot(s+1)\qquad\textrm{
for any }\qquad s\in[n]\ ,$ (127)

with $\eta_{s}(x)$ defined in Eq. (20). To obtain Eq. (126), we used that

$\displaystyle c_{j}|x\rangle$ $\displaystyle=e^{i\pi\beta_{j}(x)}\ket{x\oplus
e_{j}}\qquad\text{ for all }\quad j\in[2n]\quad\textrm{ and }\quad
x\in\\{0,1\\}^{n}$ (128)

because

$\displaystyle c_{2j-1}|x\rangle=(-1)^{\eta_{j}(x)}|x\oplus
e_{j}\rangle\qquad\text{and}\qquad
c_{2j}|x\rangle=-i(-1)^{\eta_{j}(x)+x_{j}}|x\oplus e_{j}\rangle\ .$ (129)

Eq. (126) motivates the following case distinction:

    1. (a)

$\cos^{2}(\vartheta/2)\geq 1/2$ (see Lines 8–9): Here $\ket{x}$ has higher
amplitude than $\ket{x\oplus e_{j}\oplus e_{k}}$ in the state
$U_{j,k}(\vartheta)\ket{x}$. The algorithm picks $z=x$ (Line 8) and sets
$s=\cos(\vartheta/2)$ (line 9). In particular, comparing with (126), it
follows immediately that the claims (123) and (124) are satisfied.

    2. (b)

$\cos^{2}(\vartheta/2)<1/2$ (see Lines 11–13): In this case the algorithm
ensures that

$\displaystyle z$ $\displaystyle=x\oplus e_{j}\oplus e_{k}$ by Line 11 (130)
$\displaystyle\beta$ $\displaystyle=\beta_{j}(x)+\beta_{k}(x)$ by Line 12
(131) $\displaystyle s$
$\displaystyle=e^{i\pi(\beta_{j}(x)+\beta_{k}(x))}\sin(\vartheta/2)\qquad$
$\displaystyle\text{ by Lines~{}\ref{alg:betakgeqj}
and~{}\ref{alg:sIsPhaseBeta}}\ .$ (132)

Because $\cos^{2}(\vartheta/2)+\sin^{2}(\vartheta/2)=1$ we have

$\displaystyle|s|^{2}$ $\displaystyle\geq\frac{1}{2}$ (133)

by the assumption that $\cos^{2}(\theta/2)<1/2$. To prove the claims (123) and
(124), it thus suffices to show the second claim (124). But this again follows
from (126) and the definitions (130) of $z$ and (132) of $s$, i.e., we have
$s=\langle z,U_{j,k}(\vartheta)x\rangle$.

  2. (ii)

$R=R_{j}$, $j\in[2n]$ is a reflection (see Lines 14–17): Here $R$ is
associated with the unitary evolution operator

$\displaystyle U_{j}=c_{j}\ .$ (134)

Its action on $\ket{x}$, $x\in\\{0,1\\}^{n}$, is described by Eq. (128), i.e.,
we have

$\displaystyle U_{j}\ket{x}$ $\displaystyle=e^{i\pi\beta_{j}(x)}\ket{x\oplus
e_{j}}\ .$ (135)

This state is proportional to $\ket{x\oplus e_{j}}$, showing that the choice

$\displaystyle z$ $\displaystyle=x\oplus e_{j}\qquad$ (Line 15) (136)
$\displaystyle\beta$ $\displaystyle=\beta_{j}(x)$
$\displaystyle(\text{Line~{}\ref{alg:RjBeta}})$ (137) $\displaystyle s$
$\displaystyle=e^{i\pi\beta}$ (Line 17) (138)

indeed ensures that the claims (123) and (124) are satisfied.

∎

Equipped with Claim 3.9, we can show that the algorithm $\mathsf{evolve}$ has
the desired functionality. The matrix $\Gamma_{0}$ computed in Line 4 is the
covariance matrix of the evolved state $U_{R}\Psi(d)$, whereas
$\Gamma_{1},\Gamma_{2}$ computed in Line 19 are the covariance matrices of
$U_{R}\ket{x}$ and $\ket{y}$, respectively. Thus $\mathsf{overlaptriple}$ in
Line 22 is invoked on the triple of states

$\displaystyle(\Phi_{0},\Phi_{1},\Phi_{2})$
$\displaystyle=(U_{R}\Psi(d),U_{R}\ket{x},\ket{y})\ .$ (139)

To check that the requirements of $\mathsf{overlaptriple}$ are satisfied,
first observe that

$\displaystyle u$ $\displaystyle=\overline{r}$ by Line 20 (140)
$\displaystyle=\langle\Psi(d),x\rangle$ by definition of $r$ (141)
$\displaystyle=\langle U_{R}\Psi(d),U_{R}x\rangle\qquad$ by unitarity of
$U_{R}$ (142) $\displaystyle=\langle\Phi_{0},\Phi_{1}\rangle\ .$ (143)

Furthermore, this is non-zero because $r$ (part of the input) is non-zero by
definition of the description $d=(\Gamma_{0},x,r)$ of $\Psi(d)$.

By the defining property of the subroutine $\mathsf{relatebasiselements}$,
Line 18 of the algorithm computes
$(\alpha,\gamma)\in\\{0,1\\}^{2n}\times[0,2\pi)$ such that

$\displaystyle c(\alpha)\ket{y}=e^{i\gamma}\ket{z}\ .$ (144)

We also have

$\displaystyle v$ $\displaystyle=e^{i\gamma}\overline{s}$ by line 20 (145)
$\displaystyle=e^{i\gamma}\overline{\langle z,U_{R}x\rangle}$ by (124) (146)
$\displaystyle=\langle U_{R}x,c(\alpha)y\rangle\qquad$ by (144) (147)
$\displaystyle=\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle\ .$ (148)

Because $|s|^{2}\geq 1/2$ (see Claim 3.9), this is non-zero.

We conclude from the the properties of $\mathsf{overlaptriple}$ that

$\displaystyle w$
$\displaystyle=\langle\Phi_{2},\Phi_{0}\rangle\qquad\qquad\textrm{ see
Line~{}\ref{alg:overlapEndFinal}}$ (149) $\displaystyle=\langle
y,U_{R}\Psi(d)\rangle\ .$ (150)

By construction of $y$ using $\mathsf{findsupport}$, we have

$\displaystyle|w|^{2}$ $\displaystyle\geq 2^{-n}\ ,$ (151)

see Eq. (122). In particular, we conclude that the triple $(\Gamma_{0},y,w)$
is a valid description of $U_{R}\Psi(d)$ with the desired property (121). This
is what the algorithm returns.

The runtime of the algorithm $\mathsf{evolve}$ is dominated by the runtime
$O(n^{3})$ of the algorithm $\mathsf{overlaptriple}$.

∎

We give pseudocode for the algorithm $\mathsf{measureprob}$ in Fig. 11.

1:$d=(\Gamma,x,r)\in\mathsf{Desc}_{n}$ 2:$j\in[n]$ 3:$s\in\\{0,1\\}$
4:function $\mathsf{measureprob}$($d,j,s$)  5: return
$\frac{1}{2}(1+(-1)^{s}\Gamma_{2j-1,2j})$ Figure 11: The subroutine
$\mathsf{measureprob}$ takes as input a description
$d=(\Gamma,x,r)\in\mathsf{Desc}_{n}$ of a Gaussian state $\Psi(d)$, an integer
$j\in[n]$ and a bit $s\in\\{0,1\\}$. It outputs the probability of obtaining
the measurement outcome $s$ when measuring the occupation number operator
$a_{j}^{\dagger}a_{j}$. The outcome probability does not depends on the global
phase of $\Psi(d)$ (which is determined by its reference state $x$ and the
overlap $r$), but only on its covariance matrix $\Gamma$.

######  Lemma 3.10.

The algorithm
$\mathsf{measureprob}:\mathsf{Desc}_{n}\times[n]\times\\{0,1\\}\rightarrow\mathbb{R}$
given in Fig. 11 runs in time $O(1)$. It satisfies

$\displaystyle\mathsf{measureprob}(d,j,s)=\langle\Psi(d),\Pi_{j}(s)\Psi(d)\rangle\qquad\textrm{
for all }\qquad d\in\mathsf{Desc}_{n},j\in[n],s\in\\{0,1\\}\ ,$ (152)

where $\Pi_{j}(s)=\frac{1}{2}(I+(-1)^{s}ic_{2j-1}c_{2j})$ is the projection
onto the eigenvalue-$s$ eigenspace of $a_{j}^{\dagger}a_{j}$.

###### Proof.

We denote the input to $\mathsf{measureprob}$ by $(d,j,s)$ where
$d=(\Gamma,x,r)\in\mathsf{Desc}_{n}$ is a description of a state $\Psi(d)$,
$j\in[n]$ and $s\in\\{0,1\\}$. Given the state $\Psi(d)$, the probability of
obtaining measurement outcome $s$ when measuring the occupation number
operator $a_{j}^{\dagger}a_{j}$ is given by Eq. (57). This is the output of
the algorithm in line 5 and gives the claim. Computing line 5 requires a
constant number of arithmetic operations, giving the runtime $O(1)$. ∎

We give pseudocode for the algorithm $\mathsf{postmeasure}$ in Fig. 12 and we
illustrate it in Fig. 13.

1:$d=(\Gamma,x,r)\in\mathsf{Desc}_{n}$ 2:$j\in[n]$ 3:$s\in\\{0,1\\}$
4:$p=\|\Pi_{j}(s)\Psi(d)\|^{2}\in[0,1]$ $\triangleright$ probability of
outcome $s$ when measuring $a_{j}^{\dagger}a_{j}$ 5:function
$\mathsf{postmeasure}$($d,j,s,p$)  6: $\Gamma^{\prime}\leftarrow
0\in\mathsf{Mat}_{2n\times 2n}(\mathbb{R})$ $\triangleright$ compute
covariance matrix of post-measurement state $\Psi^{\prime}$ 7:
$\Gamma^{\prime}_{2j,2j-1}\leftarrow(-1)^{s}$ 8: for $\ell\leftarrow 1$ to
$n-1$ do 9: for $k\leftarrow\ell+1$ to $n$ do 10: if $(k,\ell)\neq(2j,2j-1)$
then 11:
$\Gamma^{\prime}_{k,\ell}\leftarrow\Gamma_{k,\ell}+\frac{(-1)^{s}}{2p}(\Gamma_{2j-1,\ell}\Gamma_{2j,k}-\Gamma_{2j-1,k}\Gamma_{2j,\ell})$
12: $\Gamma^{\prime}\leftarrow\Gamma^{\prime}-(\Gamma^{\prime})^{T}$ 13:
$y\leftarrow\mathsf{findsupport}(\Gamma^{\prime})$ $\triangleright$ find $y$
such that $|\langle y,\Psi^{\prime}\rangle|^{2}\geq 2^{-n}$ 14:
$(\alpha,\vartheta)\leftarrow\mathsf{relatebasiselements}(y,x)$
$\triangleright$ $(\alpha,\vartheta)$ are such that
$c(\alpha)\ket{y}=e^{i\vartheta}\ket{x}$ 15: $\Gamma_{0}\leftarrow\Gamma$,
$\Gamma_{1}\leftarrow\Gamma(|x\rangle)$,
$\Gamma_{2}\leftarrow\Gamma(|y\rangle)$ $\triangleright$ covariance matrices
of $\Psi(d)$, $\ket{x}$ and $\ket{y}$ 16: $u\leftarrow\overline{r}$
$\triangleright$ $u=\langle\Psi(d),x\rangle$ 17: $v\leftarrow e^{i\vartheta}$
$\triangleright$ $v=\langle x,c(\alpha)y\rangle$ 18:
$w\leftarrow\mathsf{overlaptriple}(\Gamma_{0},\Gamma_{1},\Gamma_{2},\alpha,u,v)~{}$
$\triangleright$ $w=\langle y,\Psi(d)\rangle$ 19: return
$(\Gamma^{\prime},y,w/\sqrt{p})$ $\triangleright$ return a description of
$\Psi^{\prime}$ Figure 12: The algorithm $\mathsf{postmeasure}$ takes as input
a description $d\in\mathsf{Desc}_{n}$, an integer $j\in[n]$, a bit
$s\in\\{0,1\\}$ and a real number $p\in[0,1]$. Assuming
$p=\|\Pi_{j}(s)\Psi(d)\|^{2}$, the algorithm outputs a description of the
post-measurement state
$\Psi^{\prime}=(\Pi_{j}(s)\Psi(d))/\|\Pi_{j}(s)\Psi(d)\|$ when measuring the
number operator $a_{j}^{\dagger}a_{j}$ and obtaining the outcome $s$ where
$s\in\\{0,1\\}$. Here, $\Pi_{j}(s)=(I+(-1)^{s}ic_{2j-1}c_{2j})/2$ is the
projection onto the eigenvalue $s$ eigenspace of $a_{j}^{\dagger}a_{j}$.

(a) The algorithm $\mathsf{postmeasure}$ is given the covariance matrix
$\Gamma$ of a Gaussian state $\Psi$, $x\in\\{0,1\\}^{n}$ as well as the value
$r=\langle x,\Psi\rangle$ (which is non-zero). It is additionally given the
probability $p=\|\Pi_{j}(s)\Psi\|^{2}$ of observing the outcome $s$.

(b) After computing the covariance matrix $\Gamma^{\prime}$ of the post-
measurement state $\Psi^{\prime}$, the algorithm uses the subroutine
$\mathsf{findsupport}$ to find an element $y\in\\{0,1\\}^{n}$ such that
$\langle y,\Psi^{\prime}\rangle\neq 0$. (The value of this inner product is
not computed at this point.)

(c) The subroutine $\mathsf{relatebasiselements}$ is then used to find
$(\alpha,\vartheta)$ such that $c(\alpha)\ket{y}=e^{i\vartheta}\ket{x}$.

(d) The algorithm $\mathsf{overlaptriple}$ is used to compute the overlap
$w=\langle y,\Psi\rangle$.

(e) As argued in the proof of Lemma 3.11, the inner product $\langle
y,\Psi^{\prime}\rangle$ can be computed from $w$ and the probability $p$: it
is equal to $\langle y,\Psi^{\prime}\rangle=w/\sqrt{p}$. Thus
$(\Gamma^{\prime},y,w,\sqrt{p})$ is a description of $\Psi^{\prime}$.

Figure 13: An illustration of the algorithm $\mathsf{postmeasure}$, which
computes a description of the post-measurement state
$\Psi^{\prime}=\frac{\Pi_{j}(s)\Psi}{\|\Pi_{j}(s)\Psi\|}$ given a description
of $\Psi$.

######  Lemma 3.11.

The algorithm
$\mathsf{postmeasure}:\mathsf{Desc}_{n}\times[n]\times\\{0,1\\}\times[0,1]\rightarrow\mathsf{Desc}_{n}$
given in Fig. 12 runs in time $O(n^{3})$. Let $d\in\mathsf{Desc}_{n}$,
$j\in[n]$ and $s\in\\{0,1\\}$ be arbitrary. Let
$\Pi_{j}(s)=\frac{1}{2}(I+(-1)^{s}ic_{2j-1}c_{2j})$ be the projection onto the
eigenvalue-$s$ eigenspace of $a_{j}^{\dagger}a_{j}$ and let
$p=\|\Pi_{j}(s)\Psi(d)\|^{2}$. Then

$\displaystyle\Psi(\mathsf{postmeasure}(d,j,s,p)$
$\displaystyle=\frac{\Pi_{j}(s)\Psi(d)}{\|\Pi_{j}(s)\Psi(d)\|}\ ,$ (153)

that is, $\mathsf{postmeasure}$ computes a description of the post-measurement
state when measuring $a_{j}^{\dagger}a_{j}$ and obtaining outcome $s$.
Denoting the output of the algorithm by

$\displaystyle
d^{\prime}=(\Gamma^{\prime},x^{\prime},r^{\prime})=\mathsf{postmeasure}(d,j,s,p)\
,$ (154)

we further have

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

###### Proof.

We denote the input to $\mathsf{postmeasure}$ by $(d,j,s,p)$, where
$d\in\mathsf{Desc}_{n}$, $j\in[n]$, $s\in\\{0,1\\}$ and $p\in[0,1]$. For
brevity, let us denote the post-measurement state when measuring the
observable $a_{j}^{\dagger}a_{j}$ and obtaining outcome $s$ by

$\displaystyle\Psi^{\prime}$
$\displaystyle=\frac{\Pi_{j}(s)\Psi(d)}{\|\Pi_{j}(s)\Psi(d)\|}\ .$ (156)

In lines 6–11, the algorithm $\mathsf{postmeasure}$ computes the covariance
matrix $\Gamma^{\prime}$ of $\Psi^{\prime}$ according to Eq. (63). In line 13
the algorithm uses $\mathsf{findsupport}$ to find $y\in\\{0,1\\}^{n}$ such
that

$\displaystyle\left|\langle y,\Psi^{\prime}\rangle\right|^{2}$
$\displaystyle\geq 2^{-n}\ .$ (157)

Line 14 provides $(\alpha,\vartheta)\in\\{0,1\\}^{2n}\times\mathbb{R}$ such
that

$\displaystyle c(\alpha)\ket{y}=e^{i\vartheta}\ket{x}\ .$ (158)

Line 15 of the algorithm sets the matrices
$(\Gamma_{0},\Gamma_{1},\Gamma_{2})$ equal to the covariance matrices of the
three states

$\displaystyle(\Phi_{0},\Phi_{1},\Phi_{2})=(\Psi(d),\ket{x},\ket{y})\ .$ (159)

We check that the conditions for applying $\mathsf{overlaptriple}$ in Line
(18) are satisfied. We have

$\displaystyle u$ $\displaystyle=\overline{r}\qquad\qquad\qquad\textrm{ by
Line~{}\ref{alg:measureU}}$ (160)
$\displaystyle=\langle\Psi(d),x\rangle\qquad\ \ \textrm{ by definition
of~{}$r$}$ (161) $\displaystyle=\langle\Phi_{0},\Phi_{1}\rangle\ ,$ (162)

and this is non-zero because $d=(\Gamma,x,r)$ is a valid description (hence
$r\neq 0$). Similarly as before, we also have

$\displaystyle v$ $\displaystyle=e^{i\vartheta}\qquad\qquad\qquad\textrm{ by
Line~{}\ref{alg:measureV}}$ (163) $\displaystyle=\langle
x,c(\alpha)y\rangle\qquad\ \ \ \textrm{ by
Eq.~{}\eqref{eq:alphaytranslationproperty}}$ (164)
$\displaystyle=\langle\Phi_{1},c(\alpha)\Phi_{2}\rangle\ .$ (165)

In particular, this is also non-zero. The requirements to run
$\mathsf{overlaptriple}$ in Line 18 are therefore met, and Line 18 returns

$\displaystyle w$ $\displaystyle=\langle\Phi_{2},\Phi_{0}\rangle$ (166)
$\displaystyle=\langle y,\Psi(d)\rangle\ .$ (167)

It remains to show that $(\Gamma^{\prime},y,w/\sqrt{p})$ (the expression
returned by the algorithm) is a valid description of the post-measurement
state $\Psi^{\prime}$, and to establish the bound

$\displaystyle|w/\sqrt{p}|^{2}$ $\displaystyle\geq 2^{-n}$ (168)

in order to prove Eq. (155).

Inserting $\Psi^{\prime}=\Pi_{j}(s)\Psi/\sqrt{p}$, Eq. (157) implies that

$\displaystyle 2^{-n}$ $\displaystyle\leq|\langle
y,\Psi^{\prime}\rangle|^{2}=\frac{1}{p}|\langle
y,\Pi_{j}(s)\Psi\rangle|^{2}=\frac{1}{p}|\langle\Pi_{j}(s)y,\Psi\rangle|^{2}\leq\frac{1}{p}\|\Pi_{j}(s)y\|^{2}\cdot\|\Psi\|^{2}\
$ (169)

because $\Pi_{j}(s)$ is self-adjoint and with the Cauchy-Schwarz inequality.
In particular, we have $\Pi_{j}(s)y\neq 0$ and thus

$\displaystyle\Pi_{j}(s)y$ $\displaystyle=y$ (170)

since any number state $\ket{y}$ is an eigenvector of the projection
$\Pi_{j}(s)$. Inserting (170) into (169) and using (167) we obtain the bound

$\displaystyle 2^{-n}\leq\frac{1}{p}|\langle
y,\Psi\rangle|^{2}=\frac{|w|^{2}}{p}\ ,$ (171)

establishing (168). Eq. (170) and the self-adjointness of $\Pi_{j}(s)$ also
imply that

$\displaystyle\langle y,\Psi^{\prime}\rangle$
$\displaystyle=\frac{1}{\sqrt{p}}\langle
y,\Pi_{j}(s)\Psi\rangle=\frac{1}{\sqrt{p}}\langle
y,\Psi\rangle=\frac{w}{\sqrt{p}}\ .$ (172)

Since $\Gamma^{\prime}$ is the covariance matrix of $\Psi^{\prime}$ and
$p=\|\Pi_{j}(s)\Psi(d)\|^{2}$ is the probability of obtaining outcome $s$ when
measuring $a_{j}^{\dagger}a_{j}$, this shows that
$(\Gamma^{\prime},y,w/\sqrt{p})$ is a valid description of $\Psi^{\prime}$ as
claimed.

The complexity of the algorithm is dominated by $\mathsf{overlaptriple}$,
which takes time $O(n^{3})$.

∎

###  3.3 Initial states for computation

Using the algorithm $\mathsf{evolve}$, it is straightforward to generate a
description of a state that is obtained by applying a sequence of Gaussian
unitaries (generators) to the vacuum state. This is all that is typically
needed to describe initial states.

In cases where we do not need to fix the overall phase, we can generate a
description from the covariance matrix. The algorithm $\mathsf{describe}$
takes as input the covariance matrix $\Gamma$ of a Gaussian state $\Psi$ and
outputs a description $d\in\mathsf{Desc}_{n}$ of a Gaussian state which is
equal to $\Psi$ up to a global phase. It is given in Fig. 14 and it simply
uses the subroutine $\mathsf{findsupport}$ and Eq. (45).

1:$\Gamma\in\mathsf{Mat}_{n\times n}(\mathbb{R})$ covariance matrix of a pure
Gaussian state $\Psi$ 2:function $\mathsf{describe}$($\Gamma$)  3:
$y\in\\{0,1\\}^{n}$ 4: $y\leftarrow\mathsf{findsupport}(\Gamma)$
$\triangleright$ find $y$ such that $|\langle y,\Psi\rangle|^{2}\geq 2^{-n}$
5: $\sigma\leftarrow\mathsf{Pf}(\Gamma)$ 6: $r\leftarrow\sqrt{\sigma
2^{-n}\mathsf{Pf}(\Gamma(\ket{y})+\Gamma)}$ $\triangleright$ compute $|\langle
y,\Psi\rangle|$ 7: return $(\Gamma,y,r)$ Figure 14: The algorithm
$\mathsf{describe}$: Given the covariance matrix $\Gamma$ of a Gaussian state
$\ket{\Psi}$, it outputs $d\in\mathsf{Desc}_{n}$ such that