The Power of Noisy Fermionic Quantum Computation

We consider the realization of universal quantum computation through braiding of Majorana fermions supplemented by unprotected preparation of noisy ancillae. It has been shown by Bravyi [Phys. Rev. A 73, 042313 (2006)] that under the assumption of perfect braiding operations, universal quantum computation is possible if the noise rate on a particular 4-fermion ancilla is below 40%. We show that beyond a noise rate of 89% on this ancilla the quantum computation can be efficiently simulated classically: we explicitly show that the noisy ancilla is a convex mixture of Gaussian fermionic states in this region, while for noise rates below 53% we prove that the state is not a mixture of Gaussian states. These results were obtained by generalizing concepts in entanglement theory to the setting of Gaussian states and their convex mixtures. In particular we develop a complete set of criteria, namely the existence of a Gaussian-symmetric extension, which determine whether a state is a convex mixture of Gaussian states.


Introduction
One interesting route towards universal quantum computation is through the realization of Majorana fermion qubits. Such Majorana fermion qubits, encoded in pairs of non-local fermionic zero-energy modes, are believed to be present in various quantum systems such as a ν = 5/2 fractional quantum Hall system, p x + i p y superconductors and recently proposed topological insulator/superconductor and semiconducting nanowires/superconductor structures (see [1] for a review). Braiding operations on such Majorana fermion qubits can implement certain topologically protected gates, namely single-qubit Clifford gates [2,3]. A system of Majorana fermions initially prepared in a fermionic Gaussian state (see definition below) and undergoing braiding operations (which are a subset of non-interacting fermion operations 4 ) can be classically efficiently simulated if it is not supplemented with additional resources. Universal quantum computation can be achieved if this supplement consists, for example, of either (i) gates which use quartic interaction between Majorana fermions [4] or (ii) a quartic parity measurement [5] or (iii) two ancillae in certain special states |a 4 and |a 8 , involving respectively 4 and 8 Majorana fermions [2]. The advantage of the last realization is that even when these ancillae are noisy, one can purify them using braiding operations. This distillation of clean ancillae states, which are then used to implement the required gates for universality, is similar to the magic-state-distillation scheme by which Clifford group operations are used to distil almost noise-free single qubit π/8 ancillae from noisy ones [6]. For the model of computation based on braiding of Majorana fermions, distilled |a 4 states give the ability to perform single qubit π/8 gates, while distilled |a 8 states allow for the implementation of CNOT gates [2]. The noise threshold for such distillation schemes, which assume that the distilling gates and operations are noise-free, is of the order of tens of per cent, which makes them attractive. But one can also ask the converse question: how noisy are these ancillae allowed to be before one can efficiently classically simulate the entire quantum computation? This question has been addressed in the case of noisy single-qubit magic states and noisy single-qubit gates [7][8][9].
In this paper, we consider a similar question for computation using Majorana fermions. It is not hard to show (see section 5) that if the noisy ancillae are convex mixtures of Gaussian fermionic states and the computation involves only non-interacting fermionic operations, one can still classically efficiently simulate such a quantum computation. Hence we set out to develop a criterion that determines whether a state is a convex mixture of Gaussian fermionic states in section 3. We find such a criterion in the form of a hierarchy of semidefinite programs, similarly as for separable states [10]. Unfortunately, the computational effort for implementing this general criterion is too large to give decisive information for our problem of interest and we have recourse to an analytical approach for the noisy |a 8 ancilla in section 4. Even though our general criterion is not immediately useful for the problem at hand, its generality and similarity to separability criteria makes it interesting in itself.
Our work is different from previous work on criteria which determine whether a fermionic state with a fixed number of fermions has a single Slater determinant or whether a fermionic state can be written as a convex combination of states with a single Slater determinant [11,12]. The important distinction is that we fix only the parity of the fermionic state and not the number of fermions. Our goal is to extend the class of Gaussian fermionic states in a natural way, by considering states which can be written as convex combinations of Gaussian states. We call such states convex-Gaussian or 'having a Gaussian decomposition'. Even though physical states have a fixed number of fermions, Gaussian fermionic states are important approximations to physically non-trivial states, such as the superconducting BCS state. Our criterion thus intends to separate Gaussian states, and convex combinations thereof, from fermionic states with a richer structure which cannot be simply viewed as states in which fermions are paired 5 . Note that the question of having a Gaussian decomposition is also different from the question of pairing that is analyzed in [14]: Gaussian fermionic states can be paired while single Slater determinant states cannot. We hope that our results of separating convex-Gaussian states from fermionic states with a richer structure may lead to tools for understanding ground states of interacting fermion systems beyond mean-field theory.

Preliminaries: fermionic Gaussian states
We consider a system of m fermionic modes, with the corresponding creation (a † k ) and annihilation (a k ) operators (k = 1, . . . , m), respecting the Fermi-Dirac anti-commutation rules, {a j , a k } = 0 and {a j , a † k } = δ jk I . For systems in which fermionic parity is conserved, it is more convenient to use the 2m Majorana fermion operators defined as c 2k−1 = a k + a † k and are Hermitian, traceless and form a Clifford algebra C 2m with {c j , c k } = 2δ jk I . 6 Any Hermitian operator X ∈ C 2m that can be written as the linear combination of products of an even number of Majorana operators is called an even operator, i.e.
where the coefficients α 0 and α a 1 ,a 2 ,...,a 2k are real, for (c a 1 c a 2 · · · c a 2k ) † = (−1) k c a 1 c a 2 · · · c a 2k . The parity of the number of fermions is conserved by the action of an even operator X as X commutes with the fermionic number-parity operator . Thus the projector onto a pure state |ψ with fixed parity is an even fermionic operator, while |ψ has eigenvalues C all = ±1 depending on whether the parity of the number of fermions in |ψ is odd or even.
Given an even state ∈ C 2m , the correlation matrix M is a 2m × 2m real, anti-symmetric matrix with elements Real, anti-symmetric matrices such as M can be brought to block-diagonal form by a real orthogonal transformation R ∈ S O(2m), i.e.
Let us now define fermionic Gaussian states. Fermionic Gaussian states are even states of the form = K exp(−i i = j β i j c i c j ) with a real anti-symmetric matrix β i j and normalization K . Hence Gaussian fermionic states are ground states and thermal states of non-interacting fermion systems. One can block-diagonalize β and re-express in standard form as wherec = R T c with R block-diagonalizing the matrix β i j . The coefficients λ j (which can be expressed in terms of the eigenvalues of β i j ) will lie in the interval [−1, 1]. For Gaussian pure states, λ j ∈ {−1, 1} so that M T M = I , while for Gaussian mixed states M T M < I . In the theory of Gaussian fermionic states, a special role is played by those unitary transformations which map Gaussian states onto Gaussian states. These are the transformations generated by Hamiltonians of non-interacting fermionic systems, i.e. Hamiltonians which are quadratic in Majorana fermion operators. In this paper, we will refer to these unitary transformations as fermionic linear optics (FLO) transformations, as they, similarly as for bosonic linear optics transformations, have the property that where U is an FLO transformation and R ∈ S O(2m). Hence the transformations which blockdiagonalize M (or β) are FLO transformations. Note that C all is invariant under any unitary FLO transformation U , as U C all = C all U . It is clear from the standard form, equation (4), that a fermionic state is fully determined by its correlation matrix M. The expectation of higher-order correlators for a Gaussian fermionic state are efficiently obtained by virtue of Wick's theorem where M| a 1 ,...,a 2 p is the sub-matrix of M which contains only the elements M jk with j, k = a 1 , . . . , a 2 p and Pf(.) is the Pfaffian 7 . Their efficient description (2) combined with the possibility to efficiently evaluate the expectation value of observables (6) makes fermionic Gaussian states a valuable tool for approximating ground states, thermal states or dynamically generated states using (generalized) Hartree-Fock methods of interacting fermion systems, see e.g. [13,[15][16][17]. Their concise representation together with equation (5) also allows for an efficient classical simulation of quantum computations that employ only FLO operations and are initialized with Gaussian states [18][19][20][21].

Beyond Gaussian states: Gaussian decompositions
Given all the applications of Gaussian states, it is natural to try to enlarge this set of states to form a convex set, in analogy with extending product states to separable states. A first observation in this direction is that mixed Gaussian states are straightforward convex mixtures of pure Gaussian states. This is easily seen using the standard form (4), i.e. any Gaussian state can be written as , where 0 p k 1 (by setting p k = 1/2 we obtain the state I /2 m ). However, the convex mixture of two Gaussian states is, in general, not a Gaussian state [22,23]. This motivates the following definition: Definition 1 (Convex-Gaussian). An even density matrix ∈ C 2m is convex-Gaussian iff it can be written as a convex combination of pure Gaussian states in C 2m , i.e. = i p i σ i where σ i are pure Gaussian states, p i 0 and i p i = 1.

Remark 1.
Since every even state in C 2m is defined by 2 2m − 1 real coefficients, as is clear from equation (1) fixing the normalization α 0 = 1/2 m , Carathéodory's theorem guarantees that for any convex-Gaussian state there exists a Gaussian decomposition with at most 2 2m pure Gaussian states.

Remark 2.
There is both similarity and difference between the set of separable states and the set of convex-Gaussian states. For example, for separable states local unitary transformations captured by O(m) parameters relate the extreme points of the set, whereas for convex-Gaussian states the extreme points are the pure Gaussian states which are related by FLO transformations, captured by O(m 2 ) parameters. Note that the pure states σ i in decomposition of a convex-Gaussian state will typically have a correlation matrix M(σ i ) which is block-diagonal in a different {c i } basis. If we were to restrict ourselves to convex combinations of Gaussian states which are in standard form using one and the same set {c i }, we would obtain a convex set isomorphic to the set of separable states diagonal in the classical bit-string basis.
We will now prove that for systems of 1, 2 or 3 fermionic modes, all pure even states are Gaussian (a different proof of this result can be found in [24]). From this observation it is then immediate that any even state ∈ C 2m , with m = 1, 2, 3, is convex-Gaussian. In order to prove this result we start with a useful known lemma, reproduced here for completeness, which shows where S 2m is the set of all permutations of 2m symbols. The Pfaffian is non-zero only for an anti-symmetric matrix of even dimension. that any fermionic even density matrix has a correlation matrix with eigenvalues in the complex Proof. We present two proofs. The correlation matrix of an even density matrix is an antisymmetric matrix, and hence it can be block-diagonalized as in equation (3). Its eigenvalues are thus ±iλ k and λ k ∈ [−1, 1], as from equation (2) λ k is the expectation of the Hermitian operator ic 2k−1c2k which has ±1 eigenvalues. We can also prove the lemma by using a dephasing approach introduced in [25]. We provide this alternative proof as it clearly shows that M T M = I iff is a pure Gaussian state, and because we later use elements of this proof in proposition 1. Let ∈ C 2m be any density matrix and let Note that the FLO transformation U k induces an orthogonal transformation R k that leaves The iteration in k has the effect of dephasing in the eigenbasis of each ic 2k−1 c 2k . Thus m is a density matrix which contains only the mutually commuting operators ic 2k−1 c 2k , and hence its eigendecomposition involves the eigenstates of these operators which are Gaussian pure states. Note that m is not necessarily a Gaussian state, but its eigenstates are Gaussian pure states, each of which has correlation matrix block-diagonal in the same basis. This implies that the correlation matrix M is that of a Gaussian mixed state. Note that the dephasing procedure can only increase the entropy of . Therefore, if M T M = I at the end of this procedure, then the state at the end of the procedure, m , must be a pure Gaussian state. By the same token, at the beginning of the procedure (with the same M T M = I ) must be pure Gaussian. Thus M T M = I iff is a pure Gaussian state. Proposition 1. Any even pure state |ψ ψ| ∈ C 2m for m = 1, 2, 3 is Gaussian.
Proof. The statement is trivial for m = 1. Consider m = 2 and observe that when we blockdiagonalize its correlation matrix, any even pure state can written as some linear combination of I , quadratic terms ic 2k−1 c 2k and C all . The dephasing procedure in the proof of lemma 1 leaves any such even density matrix invariant-note that the operator C all is invariant under FLO transformations. As the output is always convex-Gaussian, the input state was a pure Gaussian state. The proof for m = 3 follows the same structure as for m = 2, but employs the additional observation that C all |ψ ψ| = ±|ψ ψ| for even pure states. Any pure state density matrix |ψ ψ| ∈ C 6 , when written in the basis in which its correlation matrix is block-diagonal, can be expressed as |ψ ψ| = α I + βC all + i k γ k c 2k−1 c 2k + i< j<k<l η i jkl c i c j c k c l . The condition that C all |ψ ψ| = ±|ψ ψ| implies that the quartic terms only have paired correlators c 2k−1 c 2k : this is because C all interchanges the quartic and quadratic terms. Thus |ψ ψ| has quartic terms which only involve commuting operators c 2k−1 c 2k . As argued in the proof of lemma 1, |ψ ψ| is then a convex mixture of pure Gaussian states, and the result follows.
As for separable states [26][27][28], the volume of convex-Gaussian states is non-zero for any m: any state ∈ C 2m is convex-Gaussian if it is 'close enough' to the maximally mixed state I /2 m , which is Gaussian. The following lemma proves this, using arguments similar to those in [27]: Theorem 1. For every even state ∈ C 2m there exists a sufficiently small > 0 such that Proof. Consider the set {σ ( λ, π )} of pure Gaussian states, with each σ ( λ, π) ∈ C 2m defined as where λ ∈ {−1, 1} m and π ∈ S + 2m , with S + 2m the set of permutations with positive signature (these correspond to S O(2m) rotations of the correlation matrix, i.e. FLO operations). There are (2m)!/m! pure states in this set and we will show that they form an over-complete basis for C 2m . We prove this by showing that any Hermitian can be written as where w ( λ,π) 0 and c 0. Furthermore, we can write If we find the decomposition in equation (9) with a certain c (which one could try to minimize in our construction of the decomposition), it then follows immediately that 1 1+c + c 1+c I 2 m is convex-Gaussian. Taking = 1 1+c we obtain our result. One could in principle upper bound c as a function of m. Equation (9) can be proved by considering how to express a single correlator α a 1 ,...,a 2k c a 1 c a 2 · · · c a 2k in terms of Gaussian states and I . Let us call the subset {a 1 , . . . , a 2k } = S and the coefficient α a 1 ,...,a 2k = α. Let C S = c a 1 c a 2 . . . c a 2k . We will use Gaussian mixed states ξ( λ, π) = 1 2 m m k=1 (I + iλ k c π(2k−1) c π(2k) ) where λ i ∈ [−1, 1]; such states can in turn be expressed as convex mixtures of pure states σ ( λ, π) for λ i ∈ {−1, 1} as in equation (8). We have By choosing λ k ∈ {−1, 1, 0} we can make sure that (i) |α|ξ( λ, π) = w λ,π ξ( λ, π) has the correct expectation α for the correlator C S (note that the sign of λ k can be chosen to fit the sign of α) and (ii) ξ( λ, π) has no higher-order correlators (by choosing some λ k = 0). The state ξ( λ, π) will then also have lower-order non-zero correlators C S for S ⊆ S, including S = ∅ corresponding to I . Thus we repeat the procedure with these various lower-order correlators, i.e. we represent them as a Gaussian mixed state ξ( λ, π) plus additional correlators corresponding to subsets of S . Proceeding in this way, we remove all correlators except I and end up expressing αC S = λ,π w ( λ,π) ξ( λ, π ) − β I for β = 2 −m λ,π w λ,π . This procedure is probably far from optimal, i.e. gives rise to a large c. One can optimize this by applying the procedure to a density matrix from which we repeatedly take out Gaussian mixed states, i.e. → − wξ , which remove the current highest-order correlator and the Gaussian state ξ is chosen to be the one with minimum (but non-negative) weight Tr wξ . Now that we have established the essential features of the set of convex-Gaussian states, the natural question that follows is how to determine whether a state has a decomposition in terms of Gaussian pure states. As for separability, there are two sides to this question. One is to find a Gaussian decomposition, a problem which is likely, as in the entanglement case, to be computationally hard in general. The reverse question is to find a criterion which establishes that a state is not convex-Gaussian. Note that the goal is to find a criterion which acts linearly on the pure Gaussian states, so that it can naturally be extended to convex mixtures thereof.
, will be useful in this context as it has been shown to lead to a necessary and sufficient criterion for a state to be Gaussian: Lemma 2 (Bravyi [19]). An even state ∈ C 2m is Gaussian iff [ , ⊗ ] = 0.
The operator has the important property of being invariant under U ⊗ U where U is any FLO transformation. The action of on ⊗ can be appreciated by considering the trace norm . 1 of the positive operator ( ⊗ ) : Tr where M is the correlation matrix of the state . Corollary 1 says that the state of two copies of pure Gaussian states is contained in the null space of the operator . In lemma A.1, we will prove that the null space of the operator is contained in the symmetric subspace Sym 2 (C 2 m ) spanned by vectors |ψ which are invariant under the SWAP operator P, i.e. |ψ = P|ψ . In particular, the null space of is spanned by states |ψ, ψ where ψ is Gaussian, while the symmetric subspace is spanned by |φ, φ for any φ. We will call the null-space of the Gaussian-symmetric subspace. As is a sum of mutually commuting Hermitian operators c i ⊗ c i with µ i = ±1 eigenvalues, the projectors on the eigenstates are P µ = 1 Hence the Gaussian-symmetric subspace is spanned by 2m m orthogonal eigenvectors P µ with the property that i µ i = 0. Clearly, the dimension of the Gaussian-symmetric subspace 2m m is smaller than the dimension of the symmetric subspace Sym 2 (C 2 m ), which is 2 m +1 2 for any m 1. In [10,29], the authors established a series of tests for separability based on the existence of a symmetric extension of any separable state. The usefulness of these tests relies on the fact that they correspond to semi-definite programs, which for small numbers of extensions can be implemented numerically. Furthermore, these tests are complete, in the sense that an entangled state does not have a symmetric extension to an arbitrary number of parties in the extension. Here we formulate a similar series of extension tests which are all passed by convex-Gaussian states, but non-convex-Gaussian states fail in some of them. Our criterion is based on the following simple observation: let a density matrix ∈ C 2m be convex-Gaussian, i.e. we can write = i p i σ i with σ i pure Gaussian states, then there exists a symmetric extension ext ∈ C ⊗n 2m , namely ext = i p i σ ⊗n i , which is annihilated by acting on any pair of tensorfactors, and Tr 2,...,n−1 ext = . This immediately leads to the following feasibility semi-definite program: Program 1. Input: ∈ C 2m and an integer n 2 Body: Is there a ext ∈ C ⊗n As there exists an isomorphism between C ⊗n 2m and C 2m,n (see the explicit map in [19] for C 2m ⊗ C 2m and C 4m ), one could alternatively express program 1 as an extension of to a physical system with 2m, n Majorana fermions.
Note that we do not need to enforce any permutation or Bose symmetry on the extension 8 because of the following. By lemma A.1 in the appendix, the null space of k,l is spanned by vectors |ψ, ψ k,l where |ψ is any pure Gaussian state. Thus the intersection of null spaces of all k,l is spanned by vectors |ψ ⊗n where ψ is a pure Gaussian state, which are vectors in the symmetric subspace of n parties. As ext lies in this null space, it will thus be Bose symmetric [30], i.e. π ext = ext , where π is any permutation in S n .
It is easy to see that any state has a Bose-symmetric extension; hence the stronger constraint imposed by is crucial. We say that a state has a Gaussian-symmetric extension if there exists an n for which program 1 returns yes. An explicit construction of the SDP for the case n = 2 is given in appendix B, and some numerical results are discussed in the following section. Here we will prove that only convex-Gaussian states have a Gaussian-symmetric extension to an arbitrary number of spaces, thus showing that these criteria form a complete hierarchy for deciding membership in the set of convex-Gaussian states.

Theorem 2. An even state ∈ C 2m has a Gaussian-symmetric extension to an arbitrary number of parties if and only if is convex-Gaussian.
Proof. One direction is immediate, and has already been alluded to in the formulation of the SDP. Assume that is convex-Gaussian, i.e. = i p i σ i where σ i are pure Gaussian states and { p i } is a probability distribution. Then the extension ext = i p i σ ⊗n i has the property that k,l n ext = 0, where k,l acts on the tensor factors k and l by corollary 1. To prove the other direction, we will invoke the quantum de Finetti theorem [30]. Let n ext 0 be the extension of to C ⊗n 2m such that k,l n ext = 0 and Tr 2,...,n n ext = . Using lemma A.1 this shows that n ext is Bose-symmetric. Let n 1,2 be the reduced density matrix for the first two tensor factors, i.e. n 1,2 = Tr 3,...,n n ext . Then theorem II.8 in [30] tells us that with m (n) = 4×2 m n . Here τ a (n) ∈ C 2m are density matrices and γ a (n) are probabilities summing up to 1. We have a γ a (n)τ a (n) ⊗ τ a (n) For fixed m, in the limit n → ∞ the upper bounds m (n),˜ m (n) → 0, and by equation (12) n 1,2 → a γ a (n)τ a (n) ⊗ τ a (n). In this limit, equation (13) implies, using the observation in equation (11), that for each a either γ a (n) → 0 or τ a (n) converges to a Gaussian pure state. Hence for n → ∞ the state n 1,2 converges to a convex mixture of two copies of Gaussian pure states, and its reduced density matrix converges to a convex-Gaussian state.

Remark 3.
More direct proofs (by establishing a Gaussian-fermionic de Finetti theorem directly using the existence of an extension in the Gaussian-fermionic subspace) and quantitative versions of this result should be possible.
Note that convex-Gaussian states always have a symmetric extension which is separable in relation to all tensor factors. This property can be used to provide an alternative formulation of the necessary and sufficient criterion for a state to be convex-Gaussian:

Proposition 2.
A state ∈ C 2m is convex-Gaussian iff there exists an extension ext ∈ C 2m ⊗ C 2m that is a feasible solution to program 1 -input and n = 2-with the additional constraint that ext is separable between the tensor factors.
Proof. As before one direction is immediate. For the other direction, assume that there exists a state ext ∈ C 2m ⊗ C 2m which is separable, i.e. it is of the form ext = i p i τ i ⊗ ς i , with {τ i } and { i } sets of density matrices in C 2m and { p i } a probability distribution. From the constraint ext = 0, it follows that where the notation M( ) represents the correlation matrix of the state .
Using the triangle and Cauchy-Schwartz inequalities, we have From lemma 1, it follows that TrM T M = 2 m i=1 |λ i | 2 2m with equality iff M is the correlation matrix of a pure Gaussian state. In order for equality to hold in equation (15), each τ i (and i ) must thus be a pure Gaussian state and therefore = i p i τ i is convex-Gaussian.
Including this extra separability constraint in program 1 breaks its SDP structure-checking separability can in itself be cast as an SDP [29], but this would give a nested SDP structure to the program deciding if a state is convex-Gaussian. A direct semi-definite relaxation of this criterion is obtained by employing the partial-transposition test for separability [31,32], leading to the following SDP: Here T 1 is the transposition map on the first tensor factor. As for the previous program, a negative answer means that the state is not convex-Gaussian. For the converse, it is likely that the PPT criterion together with the constraint that the extension lives in the Gaussian-symmetric subspace is not sufficient to enforce separability (e.g. there exist symmetric bound-entangled states [33]); hence a positive answer is non-decisive, but we leave this as an open question.

Applications: depolarized |a 8
As mentioned in the introduction, one way to turn ν = 5/2 topological quantum computation into a universal quantum computation scheme is to assume that one has access to the auxiliary states |a 4 ∈ C 4 and |a 8 ∈ C 8 , as shown in [2]. In this section, we want to assess the amount of noise that can be added to these extra resources before they turn convex-Gaussian. From proposition 1, we know that any pure state in C 4 is Gaussian, and therefore any noisy version of |a 4 is convex-Gaussian. We thus concentrate on the state |a 8 when undergoing depolarization. We consider the state with 0 p 1 and show that p * , the noise threshold above which a 8 ( p) is convex-Gaussian, is in the interval [ 8 15 , 8 9 ], depicted in gray. In this interval we do not yet know whether a 8 is convex-Gaussian. For p < 0.40, Bravyi [2] showed that a 8 can be distilled by noiseless braiding operations to enable topological quantum computation (TQC) with Clifford group elements (for quantum universality, one also needs a sufficiently low-noise |a 4 ancilla). Our results from the SDP on the existence of a Gaussian-symmetric extension of a 8 ( p) for p 0.25 do not provide additional information beyond these analytical results. [2] has shown that a 8 ( p) can be distilled for p < 40% (corresponding to the parameter 8 defined in [2] equal to 0.38). We would like to determine the noise threshold p * such that ρ a 8 ( p p * ) is convex-Gaussian. Note that ρ a 8 is Gaussian only for p = 1. The results in this section are summarized in figure 1.
In order to give a lower bound on p * we implemented the SDP program 1 for fixed n = 2 (see appendix B) with the aid of the cvxopt library for Sage [34]. The numerical results indicate that p * 0.25, as it is always possible to find a Gaussian-symmetric extension for a 8 ( p) with n = 2 for noise rates above this value. Our implementation of the feasibility SDP takes approximately one day for each value of p in a computer with 2.3 GHz processor, and takes about 30 GB of RAM memory. The implementation of the SDP with an additional PPT condition on the extension is likely to provide stronger results, but it is more computationally intensive. Given these limitations, we provide an alternative analysis of the properties of a 8 ( p) using the notion of witnesses: Hermitian operators which have expectation values in a restricted range (e.g. non-negative for all Gaussian states) whereas expectation values on non-Gaussian states can extend beyond this range (e.g. be negative). The existence of such witness operators is guaranteed by the fact that they represent hyperplanes separating the convex-Gaussian states from all even Hermitian operators in C 2m (in addition, similarly as in [29], witnesses could be constructed based on a negative output of feasibility SDP).
Clearly, such Hermitian witness operators, say W , should include terms which are quartic or higher-weight correlators, as the expectation value of quadratic Hamiltonians only depends on the correlation matrix of a state. For a 8 ( p) the obvious choice for such a witness operator is W = |a 8 a 8 |, for which TrW a 8 ( p) = 1 − 15 p 16 . Let us understand how we can bound min ψ g TrW |ψ g ψ g | = | ψ g |a 8 | 2 where ψ g is any Gaussian pure state. We will use the fact that the correlation matrix of |a 8 is the null matrix, i.e. for all Majorana operators c i , c j , Tr c i c j |a 8 a 8 | = 0, which follows directly from equation (17). We will now prove that | ψ g |a 8 | 2 1 2 which shows that for p < 8 15 ≈ 0.53, the state a 8 ( p) is not convex-Gaussian. The following proposition proves our claimed bound and can be intuitively understood by interpreting the state |a 8 as a maximally entangled state while Gaussian states are product states: Proposition 3. For all Gaussian pure states |ψ g , we have | ψ g |a 8 | 2 1 2 .
A tighter bound than the one in proposition 3 cannot be excluded if one uses further properties of the state |a 8 -such a tighter bound would lead to a reduction of the gray area in figure 1.
To give an upper bound on p * we construct explicit Gaussian decompositions of a 8 ( p) for large values of p. To do this we consider a subset of the convex-Gaussian states that share some key properties with a 8 ( p), namely: (i) zero correlation matrix, and (ii) only a small fraction of all possible correlators has non-zero coefficients.
To exploit these symmetries we define three types of states: where M i (λ) , for i = 1, 2, 3, is the Gaussian state generated by the correlation matrix M i (λ), with λ ∈ [−1, 1] 4 , and M i (λ) assumes one of the following three forms: By construction, these states are convex-Gaussian, have null correlation matrix, and any convex combination of them has an expansion in terms of Majorana operators with non-vanishing coefficients only on the same correlators as a 8 ( p). The aim is to find for which p's it is possible to decompose a 8 ( p) as a convex sum of these types of convex-Gaussian states.
In order to describe the convex hull of these states, first note that M 2 and M 3 are related to M 1 by permutations (different choices of pairings): Given that the transformation connecting the correlation matrices is formed by the direct sum of identical permutations, this transformation is in S O(8)-one always gets the signature of the permutation squared, and therefore the corresponding determinant is always 1. Rotations in S O of the correlation matrix induce unitary transformations on the level of the states [19]: P 12 → U 2 , P 13 → U 3 . Bearing in mind this connection among the three classes of states above, we now determine the extreme points of the set of type 1 states.

Lemma 3 (Extreme points for type 1 states). Let
Then the extreme points of S 1 are given by , with x ∈ {0, 1} 4 and ¬x is the bitwise negation of x.
Proof. It follows immediately from the expansion of the state 1 (λ): Defining in a similar fashion S 2 , S 3 , φ 2 (x) and φ 3 (x), we can then define the convex hull of any convex combination of the three types of states as S = conv{S 1 , S 2 , S 3 }. More explicitly, where the summation in x is to be carried out over its binary expansion in three bits.
The problem now reduces to determining for which values of p ∈ [0, 1] the linear system , with {γ i } a probability distribution, has a solution. Simple inspection shows that such a Gaussian decomposition is always possible whenever p 8/9 ≈ 0.89, which is then an upper bound on p * .

Discussion: classical simulation of fermionic quantum computation
In contrast to bosonic linear optics, quantum computations based on FLO can be efficiently simulated by a classical computer even when augmented by measurements of fermionic number operators [18][19][20]. The simulation relies on the following facts With these three ingredients, the evolution of any initial Gaussian state can be efficiently followed by a classical computer at any time in the computation, and the output distribution can be exactly evaluated. In [19,21], the allowed transformations in (ii) were extended to noisy Gaussian maps. These are completely positive channels, not necessarily FLO unitaries, that map Gaussian states onto Gaussian states, and as such still admit efficient classical simulation.
Our results imply that if p 8/9, a 8 ( p) is convex-Gaussian and this allows for the classical simulation of the computation as follows. Given a noise strength p 8/9, one finds the decomposition of a 8 ( p) in S, by solving the linear system a 8 -this has to be done only once for a fixed p. As each φ i (x) is in itself a convex sum of two pure Gaussian states, this leads to a decomposition of a 8 ( p) in terms of at most 48 pure Gaussian states. Let p i be the probability associated with the ith pure Gaussian state in such a decomposition. Then, whenever a state a 8 ( p) is requested in the computation, one samples from { p i } and chooses the corresponding pure Gaussian state. From this point onwards, the simulation follows the scheme summarized in (i), (ii) and (iii) above.
The results in figure 1 leave some gaps in our understanding, in particular about the precise value of p * . In order to improve the lower bounds on p * one could consider distillation processes beyond the restricted set of braiding operations, but that still map Gaussian states onto Gaussian states, leaving thus the set of convex-Gaussian states unchanged. The distillation protocol of depolarized GHZ 4 states, which are isomorphic to a 8 ( p), proposed by Dür and Cirac at first sight suggests that our classical simulation threshold is tight, as it can distil a perfect GHZ 4 provided that p < 8/9. Nevertheless, their protocol uses non-FLO operations, as it assumes that one can locally create the GHZ 4 and uses distilled two-qubit maximally entangled pairs for distribution. The protocol devised by Murao et al in [35] directly distils GHZ 4 from depolarized copies of it for p < 0.7328. Despite the fact that this direct approach seems more related to our question, it employs gates that are not immediately translated into braiding operations or even to FLO transformations. Although the choice to distil via braiding operations is physically motivated-due to their topological protection-a different distillation threshold from FLO operations would suggest that either Bravyi's distillation protocol can be improved to higher depolarizing noise rates or that there exists a 'transition zone' in which the computation by braiding of Majorana fermions is neither quantumly universal nor can be efficiently simulated classically.
Note that we have explicitly kept the permutation-symmetry constraint in this program. In this appendix, we convert this problem into a feasibility semi-definite programming problem. To do that we must show that the above problem can be cast as the standard form of SDPs, namely where x ∈ R d is the unknown vector to be optimized over, c is a given vector in R d , and {F i } i=0,...,d are given symmetric matrices. For the equality constraint, A ∈ R p×d is a given matrix, with rank(A) = p, and b ∈ R p . The restriction on the rank of A demands that the linear system has at least one solution, and all the rows are linearly independent. The first step is to associate an Hermitian operator M j with each term i k c a 1 c a 2 · · · c a 2k in the expression (1). The set {M j } 0 j 2 2m−1 −1 spans all the even Hermitian operators of C 2m . Then any state ∈ C 2m can be written as Therefore, β 0, j = α j /2 m . In this way, a further 2 2m−1 coefficients are determined, and thus there remain 2 2m−3 (2 2m − 2) free parameters. Note that the normalization of ext is already guaranteed as long as Tr( ) = α 0 = 1.
The question of whether there exists an assignment of the remaining free coefficients respecting the last two constraints can be immediately posed as a feasibility SDP: minimize 0 subject to ext 0 ext = 0.
Forgetting for the moment that the correlators are in general complex, a direct correspondence with the SDP standard form can be done as follows: The equality constraint then reads This can be cast in the form Ax = b, by writing A = ([ F i ] 1 i 2 2m−3 (2 2m −2) ) and b = −[ F 0 ]. The notation [G] means a column matrix constructed by stacking each column of G below each other.
To finish the construction we must take into account that we are possibly dealing with complex values. To do that we employ the following well-known result. Proof. (⇒) Let z = x + iy be an eigenvector of Z with eigenvalue λ. Then it is easy to see that (x, y) T is an eigenvector of T(Z ) with eigenvalue λ. (⇐) By the same token, let (x, y) be an eigenvector of T(Z ) with eigenvalue λ. Then z = x + iy is an eigenvector of Z with eigenvalue λ.
The final dimension of each F i is then 2 2m+1 × 2 2m+1 , and there are 2 2m−3 (2 2m − 2) of such matrices. For the equality constraint we write The final matrix (Re(A), Im(A)) T has dimension 2 4m+1 × 2 2m−3 (2 2m − 2). Most likely, the rank of this matrix is not equal to the number of rows. Therefore, before plugging it into the SDP one must evaluate its reduced-row form (echelon form), and remove the all-zero lines.
If one wishes to include the PPT test, as in program 2, this can be done by rewriting the two positivity constraints, ext 0 and T 1 ext 0, as a single one, ext ⊕ T 1 ext 0. One thus needs to map each complex F i in the positivity constraint to F i = F i ⊕ F T 1 i . The final dimension of each real F i is then 2 2m+2 × 2 2m+2 . The equality constraint remains the same.