Some questions concerning Majorana fermions in 2D (p + ip) Fermi superﬂuids

Most of the discussion in the literature of the Majorana fermions (M.F.’s) believed to exist in so-called 2D p + ip Fermi superﬂuids, and in particular of their possible application in topological quantum computing (TQC) has analyzed the problem using the familiar Bogoliubov–de Gennes (mean-ﬁeld) formalism, which conserves total electron number only mod 2. We raise the question: Suppose that we require that electron number be conserved, period, at all stages of the calculation, then what are the consequences for the characteristic properties of M.F.’s, in particular for their (physical) braiding statistics? While we do not provide a deﬁnitive answer to this question, our considerations suggest that these consequences could be enough to destroy the possibility of using these excitations for TQC, and certainly indicate that the problem deserves further study.


Introduction
One of the most exciting ideas in the field of quantum information in recent years has been that of topologically protected quantum computation: that by burying the delicate information required, at the intermediate stages of a quantum computation, about the relative phases of the different quantum states involved in various hard-to-access degrees of freedom of a many-body system one may be able to protect all or at least part of it from the ubiquitous decoherence which is the bane of such computation. (For a general review of the idea, see the review article [1].) One particular implementation of this notion, which while not offering total topological protection is believed to provide a very useful partial degree of it, involves physical systems some of whose elementary excitations are socalled Ising anyons, that is, excitations whose properties under exchange are not the standard bosonic or fermionic * Correspondence: aleggett@illinois.edu 2 Dept. of Physics, University of Illinois at Urbana-Champaign, 1110 W. Green St, Urbana, IL 61801, USA 3 Shanghai Center for Complex Physics, Shanghai Jiaotong University, Shanghai, People's Republic of China Full list of author information is available at the end of the article ones but rather constitute a nontrivial representation of the braid group (see ref. [1,Sect. 11.A.1]).
A particular and much-studied example of Ising anyons is the Majorana fermions (hereafter MF's) which, at least within the standard mean-field (Bogoliubov-de Gennes) approach to Cooper pairing in Fermi systems, occur in a class of superconductors (or neutral Fermi superfluids), the so-called "topological superconductors (superfluids)" characterized by an effective dimensionality of 2 and a particular type of order parameter (relative wave function of the pairs, hereafter abbreviated OP. namely the so-called "p + ip" form). It seems virtually certain that both the A and B phases of superfluid liquid 3 He, provided they can be confined so as to be sufficiently 2-dimensional, would fall in this category, and until recently the majority opinion was that the bulk metallic superconductor (singlelayer) strontium ruthenate did so too (however, recent experiments have thrown doubt on this assumption). Other strongly layered superconductors such as UTe2 have also been considered as candidates. In such bulk 2D systems the theoretical prediction is that MF's may occur at points where the OP goes to zero, namely either at the boundary or towards the center of a quantum vortex. However, while there have been several claims to have produced ex-perimental evidence for them in these systems, there is at the time of writing no general agreement on this. A closely related system which in recent years has been a preferred candidate is a hybrid semiconductor-superconductor system; that is, a thin semiconducting (e.g. ln) wire with strong spin-orbit coupling, placed in a magnetic field and next to a "trivial" (e.g. BCS s-wave) superconductor. In this situation the theoretical prediction is that a type of superconducting state closely similar to the bulk p + ip one will be induced in part of the semiconducting wire, and that at the boundaries between the normal (N) and superconducting (S) regions of the wire MF's may form. In part because of the large number of control parameters which can be tuned (magnetic field magnitude and direction, spin-orbit coupling, magnitude of s-wave superconducting gap. . . ), this idea has motivated a large amount of experimental work over the last decade, starting with the 2012 paper [2] of the Delft group. Most of these papers have measured the current-voltage characteristics of the wire as a function of the various parameters and compared them with the theoretical predictions for states showing MF's. Unfortunately, much of the argument has been of the "what else could it be?" variety, and the conclusions of many of these papers are controversial. (For a rather blistering commentary, by a leading worker in the field, on the standard of evidence used in some of the papers in question, see ref. [3].) Given this state of affairs, in 2018 the (US) National Science Foundation proposed a "grand challenge" to experimental groups: not just to find a set of putative Majorana fermions, but to demonstrate explicitly that their braiding characteristics are of the Ising-anyon type; the thinking was that such a demonstration would convince even the most diehard sceptics of the validity of the concept. At the 2018 April MRS meeting several groups took up this challenge, but to the best of the present authors' knowledge none has so far claimed definitive victory. So it seems that the jury is still out. . .
In this situation it may make sense to ask whether the theoretical arguments in favor of MF's with the "standard" properties are as sound as they are usually deemed to be. Let us make it clear right away that we do not suspect any "error" in the theoretical papers, in the sense that given their starting point for the description of fermionic excitations in topological superconductors, which was, at least until very recently almost without exception the standard Bogoliubov-de Gennes (BdG) equations, we see no particular reason to suspect that the ensuing argument in any of these papers is incorrect. Our doubts refer to the BdG equations themselves (or more precisely the objects which they are taken to describe); one can ask, first, do these equations themselves need correction? secondly, if they do, do the corrected versions still produce excitations which are qualitatively of the nature of those discussed in the current literature? and thirdly, do these "revised" MF's possess all the properties claimed for the "unrevised" version, in particular do they behave under braiding as Ising anyons ? We should emphasize that we see no reason to doubt the ideas concerning Majorana fermions "in the abstract" which have been developed over the last two decades and are reviewed e.g. in ref. [1]; it is only the specific application to topological Fermi superfluids which we believe may be questionable. Thus, for example, we believe that the various calculations which have been done for spin problems such as the Kitaev honeycomb model [4] are probably entirely correct; however, the "mapping" which is sometimes claimed between these systems and the topologicalsuperfluids is actually only to the "BdG" version of the latter and thus does not evade our criticism.
Since we [5][6][7], among others (e.g. [8]) became worried by this point some years ago, and particularly in the last two or three years, a number of papers have appeared in the literature which attempt to address it. While to the best of our knowledge none of them answer the more restricted question which we shall pose in the present paper, they are certainly relevant to it, and we will comment on some of them in Appendix B.
The aim of the present paper is indeed quite modest: Given the original analysis of topological quantum computing schemes based on bulk (2D) p + ip Fermi superfluids within a mean-field (BdG) scheme, how if at all are the conclusions affected by the simplest modification of the scheme which satisfies conservation of total particle number? We emphasize that except in Appendix B we do not address either 1D systems or "measurement-based" schemes; throughout the main text of the paper, the word "braiding" refers exclusively to "physical" braiding, that is the operation of moving anyons around one another in real physical space. Whether the results obtained in the recent literature concerning 1D systems and/or measurementbased schemes have any relevance to our conclusions, or vice versa, will be discussed in Appendix B.
In Sect. 2 we briefly recap the "standard" (mean-field) approach to MF's and their physical braiding in a 2D (p + ip) Fermi superfluid; rather than treating this system as a special case of the more general theory of superconductivity, we specialize to it from the start, a policy which we hope some readers new to the specific topic may find illuminating. (We note that some of the resultant formulae, e.g. eqns. (35), may differ by minus signs etc. from the familiar textbook ones for the BCS case.) In Sect. 3 we first, in Sect. 3.1, briefly discuss an ambiguity concerning the ground state, then, making the standard assumption about the latter, introduce the simplest modification of the meanfield theory which respects overall particle number conservation; in Sect. 3.3 of this section we note the consequences of such an adaptation for various properties of the notional MF's other than braiding. Section 4 is the crux of the paper: we examine the consequences of the corrected theory for physical braiding of the MFs, with as we shall see somewhat ambiguous conclusions. Section 5 is a brief summary and conclusion.
Some of the material of this paper has been presented in various talks over the last few years, most extensively at the 2015 Capri spring school, in the thesis [5] of one of us, in an unpublished online preprint [6], and in skeleton form, in ref. [7]. Some of the content of Sects. 2.4 and 3.3 is adapted from ref. [9].

Groundstate wave function in free space
We consider a collection of spinless fermions contained in a two-dimensional (2D) volume V and subjected to a chemical potential μ which will be taken to be positive definite (cf. below) and (for the moment) to periodic boundary conditions on the many-body wave function. The total number of particles N then fluctuates, but we can define its average N and consider the thermodynamic limit V → ∞, N → ∞ N /V → const. ≡ n. Single-particle momentum eigenstates are labelled by a 2D wave vector k = p/ and created by the standard fermionic operators a + k , and the normal groundstate is a 2D Fermi sea k a + k |vac of radius k F = (4πn) 1/2 , where |vac denotes the vacuum state. We consider a BCS-like effective Hamiltonian of the form where ξ k ≡ 2 k 2 2mμ and the pairing potential has the "pwave" form where k 0 is a phenomenological upper cutoff; we assume that any particle-hole asymmetry around k F can be neglected. A straightforward generalization of the original BCS method then leads to a mean-field Hamiltonian of the form with the complex "gap" function k given by the so-called "p + ip" form 1 1 We choose the sign of the helicity to agree with ref. [10].
In the following it is crucial to bear in mind that k is assumed to have the form (4) not just close to the Fermi surface, but all the way down to k = 0 (and up to k = k 0 ), otherwise some of the consequences below cannot be drawn. The problem defined by eqns. (3) and (4) has been studied extensively in the influential paper [10], and in this subsection we mainly summarize some of the results of that reference. (We do not reproduce some of the considerations related to the connection with the "ν = 5/2" quantum Hall effect, nor those concerning the topological transition at μ = 0, since these are irrelevant in the present context.) In the BCS particle-nonconserving formalism the groundstate has the form 2 where we can choose the phases so that where ϕ k is the angle of k on the Fermi surface, so that the total mean particle number, which is twice the sum over k of |v k | 2 is just as it should be It is interesting to consider briefly the question of the "total angular momentum" of the state (5) with the choice (6). Actually, it is not clear that to the extent that we impose (physically unrealistic) periodic boundary conditions this question really makes much sense, but it is intuitively tempting to define a "pseudo-angular momentum" L ps by where ∇ k is the appropriate continuum limit of a finite difference. It is clear that the expression under the sum in the expression (7) is just |v k | 2 so L ps is exactly -N /2. In the 3D case, whether this result has any real physical meaning, and in particular whether the real physical angular momentum in an appropriate geometry, of the superfluid A phase of liquid 3-He (routinely described in bulk by the 3D version of eqn. (6)), is indeed -N /2, is a nearly 50year-old conundrum with to this day no universally agreed solution.
Since the wave function (5) does not correspond to a fixed particle number, we cannot use it directly to write down the true (particle-conserving) many-particle wave function in coordinate space. However, irrespective of the form of the u k 's and v k 's we can project BCS on to the Nparticle subspace, thereby obtaining where apart from a normalization factor with the coefficient c k given by It is interesting to rewrite the operatorĈ † in terms of the coordinate space creation operatorŝ where K(r -r ) is the Fourier transform of the coefficient c k . We expect that the long-distance behavior of K(r -r ) will be determined by the behavior of the quantity c k in the limit |k| ⇒ 0. In the simple s-wave case k = const. ≡ the latter is rather boring: c k simply tends to the constant 2μ/ giving a (rather complicated) behavior of K(r -r ) which however has no interesting long-range behavior. However, for a p-wave gap given by eqn. (4) we find that the limit of c k as k → 0 is 2μ/ 0 (k x +ik y ) so that at large values of ρ = r -r , the real-space function K(ρ) has the form K(ρ) = const.(ρ x +iρ y ) -1 (13) and the many-body wave function N (r 1 , r 2 , . . . , r N ) is apart from a constant the Pfaffian of this expression. It was the similarity of this wave function to that believed to be the form for the ν = 5/2 quantum Hall effect which led the authors of ref. [10] to suggest that not only the groundstate but also the fermionic excited states of this system might have interesting topological properties. Let us emphasize once more that eqn. (13) follows only if the form (4) of the gap function k can be extrapolated right down to k = 0. One general point which we need to note for future reference is that the coefficient c k in the many-body wave function given by eqn. (10), is simply related to the Fourier transform F k of the two-body reduced density matrix of the system (for details see e.g. ref. [11, ch. 5]:

The uniform case: spontaneous symmetry breaking and Bogoliubov quasiparticles
While there is nothing wrong with the Pfaffian (numberconserving) wave function, most textbooks on superconductivity develop the theory by relaxing number conservation and writing the groundstate in the BCS form (5): in the number basis this takes the form where N is an appropriate number eigenstate (in our case the Pfaffian) and the complex coefficients q N are peaked strongly around the mean N and related in phase (argument) by Thus the invariance of the original Hamiltonian (1) under the global U(1) symmetry transformation N → exp(iα) · N is spontaneously broken in the BCS groundstate ("spontaneously broken U(1) symmetry", hereafter abbreviated SBU(1)S).
Let us now consider the simplest many-body states containing an odd number of particles. To do so, we note that the BCS groundstate is actually a special case of a more general many-body state of the form k>0˜ k , where˜ k is now an arbitrary state vector within the 4D "occupation space" of the pair of single-particle states k, -k spanned by the basis vectors where for example |01 k is a shorthand for the state in which the single-particle state |k is empty and the singleparticle state | -k is occupied. As indicated in eqn. (5), the BCS groundstate contains for each k a particular linear combination of |00 k and |11 k ; there exists a second linear combination orthogonal to this ("excited pair" state) which likewise has even total number parity, but we shall not be interested in this in the present context. Rather, we shall be interested in the odd-number-parity ("broken pair") states |01 k and |10 k . How to build an odd-number-parity manybody state incorporating one or other of these states (let's say for definiteness |10 k )? Needless to say, there is a trivial solution: simply exclude the k in question from the product (5) and substitute the factor |10 k ! But this is not very useful, since in real life the state |10 k is likely to be generated either by thermal excitation or by some kind of external perturbation, in either case starting from the a many-body state where only the ground pair is occupied. So how do we generate |10 k from the pair groundstate k (eqn. (5))? Two obvious possibilities are However, while these are perfectly good definitions of |10 k , the operations in question are not reversible. For this and other reasons the standard procedure is the one written down by Bogoliubov and Valatin, which combines the two operations (a) and (b): and similarly with |u k | 2 + |v k | 2 = 1.
With these standard choices, the α k 's have the same fermionic commutation relations as the a k 's. The point of the above paragraph, which to anyone familiar with textbook BCS theory may seem rather trivial and superfluous, is to emphasize the point that, at least in this simple homogeneous case, it is not the fermionic excitations which are complicated but rather the "bosonic" groundstate. In particular, note that while it seems natural to look at the two terms in eqns. (19) and (20) as corresponding respectively to the creation of an extra "electron" (or "particle") and an extra "hole" respectively, so long as we start from the BCS groundstate (5) the distribution of total particle number in the two branches, while not exactly the same, differs only by a term of order n -1 fl where n fl is the rms fluctuation of total number in the groundstate, and thus vanishes in the thermodynamic limit. This point becomes more important when we go on to treat the spatially inhomogeneous case in the next subsection.
Before leaving the homogeneous case, we note for future reference two other "textbook" results: (1) the states |10 k and |01 k are certainly eigenstates of the BCS reduced Hamiltonian, since they are manifestly eigenstates of the kinetic energy and the pairing term acting on them gives zero; (2) there are two further linearly independent combinations of the a's and a † 's which one can form, namely and that each of these simply annihilates the state (5) (i.e. when acting on it, produces a state of zero norm); we will call them "pure annihilators".

Spatially inhomogeneous case: the Bogoliubov-de Gennes equations
The moment we have to deal with any situation more complicated than the simple spatially uniform one considered by BCS, the simple procedure sketched in Sect. 2.2 is liable to fail and one needs a more flexible method, which was supplied for the case of simple s-wave pairing by Bogoliubov and de Gennes soon after BCS's original work; see ref. [12, ch. 5], for a detailed account. Since then, it has been generalized to "exotic" pairings, including the p + ip case of present interest, and is now the method of choice in virtually all the existing literature on the subject. In view of the special focus of interest of this paper we shall approach the "BdG" equations which are the standard outcome of the application of this method by a rather indirect route, as follows.
Consider as previously an even number N of spinless fermions which move at T = 0 in a volume V subject to the standard periodic boundary conditions, take the thermodynamic limit and following Yang [13] assume that the groundstate has the property that there is one and only one eigenvalue of the 2-particle density matrixρ 2 which is of order N (as discussed in ref. [14], this condition is slightly weaker and more general than the property of ODLRO which Yang himself emphasizes). Then as he shows, it must always be possible to choose two sets of single-particle wave functions {m}, {m} (not necessarily eigenstates of the single-particle part of the Hamiltonian) such that we can write the eigenfunction F(r, r ) corresponding to the macroscopic eigenvalue (the "Cooper pair wave function") in the form where the relation between m andm is one-to-one. Now we saw that in the simple BCS case, for which m andm are plane-wave states k and -k, the many-body wave function was of the form (9) with (10) and (12), where the coefficient c k is related to F k by eqn. (14). Consequently, it is very tempting to assume that in the general case described by eqn. (22), while we may not know the true many-body wave function N in detail, a good "effective" form of the latter which is adequate for the description of the groundstate and the simplest elementary excitations iŝ and K(r -r ) in eqn. (12) is replaced by where in the solution for c m the negative sign in front of the square root must be taken. Thus, apart from an overall normalization constant, we can write for the number-conserving groundstate or in the BCS number-nonconserving scheme We can take the analogy with the simple BCS-like case of Sects. 2.1-2.2 a little further: There are two linear combinations of the four operators a m , a + m , am, a + m , which are pure annihilators, namely and it is immediately clear that any linear combination of α m 's and αm's will likewise annihilate the many-body groundstate. However, there is now in general an important difference with the uniform case discussed earlier: consider the state generated by applying α + m (or α + m ) to the groundstate (27). This state is manifestly an odd eigenstate of numberparity. However, in distinction to the uniform case, it is not in general an eigenfunction of the Hamiltonian or even of the mean-field approximation to it! So what we need to do is find a set of operators γ † i which are linear combinations of the α m 's and α + m 's, and which further satisfy the relation whereĤ mf is the "mean field-approximated" Hamiltonian, that is, the expression obtained by replacing the term quartic inψ(r) andψ † (r) by a sum of the Hartree and Fock terms plus terms of the form where the anomalous averages must be computed selfconsistently. If we define this procedure just yields the standard form of the BdG equations. For a Hamiltonian of the form with the above mean-field approximation made and the complex quantity (r, r ) defined by the BdG equations for the coefficients u(r), v(r) occurring in the Bogoliubov quasiparticle operator 3 are of the form (The occurrence of the integral operator in the off-diagonal term may look unfamiliar, as in much of the literature on p-wave pairing it is approximated by a gradient.) Needless to say, these equations follow directly from a mean-field decomposition of the original Hamiltonian (31), but the above derivation suggests that their validity may extend beyond this approximation. An important property of the p-wave BdG equations (35), which is slightly different from its form in the s-wave case, is that if the spinor u i v i is a solution with positive eigenvalue E i > 0, then the spinor +v * i +u * i is also formally a solution, with negative eigenvalue -E i < 0. In much of the literature the corresponding operatorγ i is said to create a negative-energy eigenstate; however, this is misleading, since it can be seen from eqn. (30) to be a linear combination of pure annihilators and thus itself just a pure annihilator.

Majorana fermions: definition and physical nature
In particle physics a "Majorana fermion", as defined in Ettore Majorana's original paper, is a fermion which is its own antiparticle. Correspondingly, in condensed matter physics an M.F. is defined as a fermionic excitation (i.e., a solution of the BdG equations) whose creation operator γ † 0 satisfies the relation In the case of a Bogoliubov quasiparticle described by eqn.
(34), this immediately implies that and by the relation mentioned in the last paragraph of Sect. 2.3, it then immediately follows that such an excitation has zero energy relative to the groundstate: However, neither (37) nor (38) nor even their conjunction, is sufficient to uniquely define an M.F.: as noted in ref. [7], both these equations are satisfied by a Bogoliubov quasiparticle in a (3D) Fermi superfluid with an order parameter (gap) proportional to say k x -k y with a k-vector exactly at a node of the gap, but we would not normally call such an excitation a MF. What we need to add is that the components of the spinor corresponding to a MF must be localized in space (with respect to at least one dimension, cf. below). What is the physical nature of a MF? [9] Intuitively, one might think that since positive-energy solutions of the BdG equations correspond to real odd-number parity manybody states, while negative-energy ones simply correspond to pure annihilators, a zero-energy solution should in some sense correspond to half of each, and this is essentially true. To see this, we first note that the special case of eqn. (29) corresponding to a MF, namely when applied to the groundstate, has not one but two possible interpretations: (1) γ † 0 generates a real odd-number-parity many-body state with zero energy relative to the groundstate, (2) γ † 0 annihilates the groundstate. We will now show that γ † 0 cannot correspond to either of these possibilities alone. To do so, let us rewrite γ † 0 , regarded as doing (1), in terms of the m-representation by: where α † m , α † m are given by the Hermitian conjugates of eqns. (28). We then see that since the operators α m , αm, which have no simultaneous solution, and thus γ † 0 cannot simply do (1). A similar argument shows that it cannot simply do (2). Hence, it must be a linear combination of the operations (1) and (2): in other words [9], a MF is a linear combination of an operator γ † 0,Bog which creates a zero-energy fermionic state (Bogoliubov quasiparticle) and a pure annihilator.
How then do we create a zero-energy Bogoliubov quasiparticle, period? It is clear that such an excitation must be a linear combination of the form (40) of the operators α + m , α + m which satisfies the BdG equations with eigenvalue zero, but as we have seen it cannot be self-conjugate, indeed it must satisfy the standard Fermi anticommutation relations. The difficulty is resolved by noting that by combining two mutually orthogonal MF solutions with an appropriate relative phase π/2 we can fulfill both conditions simultaneously: The Hermitian conjugate operator γ 0,Bog is then a pure annihilator, as it should be. The above maneuver requires that for any given system characterized by an order parameter which tolerates MF solutions, the number of such solutions should be even, thus giving rise to an integer number of zero-energy Bogoliubov quasiparticles. While the general proof that this is so requires some quite sophisticated topological considerations, let us explore it for the specific case of most direct physical interest in the present context, namely vortices in a spinless Fermi system where the Cooper pairs form in a (p + ip) state (or what is in this context equivalent, "halfquantum" vortices in a spin-1/2 Fermi superfluid with this order parameter).
For pedagogic reasons we first briefly digress to discuss the case of a simple Abrikosov(-like) vortex in a neutral BCS s-wave superfluid. Then the following argument, based on the consideration of conservation of (apparent) orbital angular momentum L relative to the vortex core, indicates that no MF solution occurs: (1) lf a subscript u denotes the angular momentum associated with the u component of the solution to the BdG equations, etc., and we assume that all angular momenta are integral, then the BdG equations imply (2) On the other hand, the self-conjugacy requirement implies that L u = -L v . (3) Hence L must be an even integer. However, for s-wave pairing the only angular momentum associated with delta is the center-of-mass term, which around the core of an Abrikosov vortex is ±1. (4) Hence we reach a contradiction. For a (p + ip) Fermi superfluid the above argument says that a MF solution is allowed, with an angular momentum L = 0 or (±)1 depending on whether the "intrinsic" angular momentum associated with the (p + ip) pairing is parallel or antiparallel to the vorticity. Moreover, an explicit calculation (e.g. [15]) then shows that for a single isolated vortex exactly one such solution indeed occurs. However, it then "miraculously" turns out that either the total number of vortices must be even, or if it is odd then a second MF must occur at the boundary of the system; thus in either case the total number of MF's must be even, giving an integer number of zero-energy Bogoliubov-fermions.
What is peculiar about these Bogoliubov fermions is that they are "delocalized", that is, their creation operators γ † 0,Bog (which from now on I shall write simply as α † 0 ) have support from two different regions of space (the neighborhood of the individual vortices) which may be widely separated. Moreover, in the case of a number of vortices greater than 2 the association of a "particular" Bogoliubov fermion with a particular pair of vortices is arbitrary: e.g. for the case of 4 vortices we could associate one Bogoliubov fermion with vortices 1 and 2 and B.F. 2 with vortices 3 and 4, but we could equally well define one of the B.F's as associated with vortices 1 and 3 and the second with 2 and 4; and so on. These features give rise to what is perhaps the most fascinating property of the MF's (or the BF's made out of them), namely their behavior under braiding, and this is the subject of the last subsection of this section.

Braiding of Majoranas
In a widely-cited 2001 paper [16], Ivanov discusses, within the standard mean-field ansatz, the braiding properties of half-quantum vortices in a 2D (p + ip) Fermi superfluid (or equivalently of vortices in a spinless superfluid), reaching the conclusion that they obey the statistics characteristic of Ising anyons. His argument essentially consists of three stages: (1) the demonstration that Majorana fermions exist on such vortices, and can be associated to form zeroenergy fermions in the standard way (2) the demonstration that the result of an interchange of two vortices is described by his eqn. (8) (see below) (3) the demonstration that when eqn. (8) is applied to a system containing more than two vortices, the resulting interchanges are described by the braid group. Since stage (1) is a relatively straightforward application of the BdG equations (see also below), and stage (3) does not seem controversial, we will concentrate here on stage (2). This needs a little discussion: First, we will be interested only in the case of interchange of two vortices (not of a vortex and antivortex, which in effect are distinguishable particles). Further, in contrast to the bulk of the literature (including ref. [16], we will explicitly consider the case in which the helicity of the vortex is opposite to that of the intrinsic angular momentum of the pairs. The reason is that in this case conservation of "local" angular momentum in the BdG equations, combined with the self-conjugacy condition, implies that the coefficients u(r) and v(r) occurring in the M.F. 's have winding number zero around the vortex core (as distinct from ±1 in the usually considered case), i.e. are independent of the angle with respect to it, which removes some irrelevant complications. Finally, we state precisely the question to which Ivanov's eqn. (8) is the answer, namely: Consider two vortices which can be either (a) in a state with no Majorana fermions (so that the number parity of the corresponding many-body state is even) or (b) in a state where each vortex carries a Majorana fermion (so that the number parity of the many-body state is odd). Imagine that we adiabatically interchange these two vortices, say clockwise; let the phase accumulated by the many-body wave function in states (a) and (b) be respectively ϕ a and ϕ b . What is the difference ϕ aϕ b of these two phases? Note that we are not interested in the "absolute" value of the phase in either state.
The derivation of eqn. (8) in the original paper [16] involves making a "cut" in the definition of the order parameter (r). The following argument, while somewhat less general, avoids the necessity of this maneuver and is essentially equivalent (more details are given in Sect. 4): Rather than considering an exchange of the two vortices, we site one of them at the center of a cylindrically symmetrical disk and consider "encircling" it with the other at some radius R which is the radius of the vortex core but the radius of the disk. The phase difference on exchange is then the square root of that acquired on encirclement. Denote the angle made by the core of the moving vortex with the origin at any given point in the encirclement process by θ 0 We assume, as Ivanov implicitly does, that the only relevant effect of the motion of θ 0 is to change the value of the complex gap delta "due to" vortex 1 at the position of vortex 2 and vice versa; we will return to this point in Sect. 4. Now, quite generally, the phase accumulated in a closed-loop adiabatic process is the sum of two terms, usually known respectively as the "monodromy phase" and the "Berry phase" (see for example ref. While each of these phases separately may depend on the choice of gauge, their sum does not. We will choose a gauge in which the phase of the "overall" OP at the moving vortex evolves "clockwise" by 2π during the encirclement owing to the effect of the stationary one and vice versa. If we consider not the absolute value of either ϕ a or ϕ b but rather their difference, the monodromy phase will be given by that accumulated by u(r, θ 0 ) and v(r, θ 0 ). Since the problem defined by the BDG equations is isomorphic to that of a spin-1/2 particle in a magnetic field which rotates through 2π in the xy-plane, with u and v corresponding to the +z-components of the spin, we find that as θ 0 moves from 0 to 2π , u and v each evolve oppositely through π , i.e. we have (Note that this dependence of u and v on θ 0 should not be confused with that on the angle with respect to their "own" vortex core, which remains constant throughout.) So the upshot is that the monodromy phase is π . As to the Berry phase, the condition u i (r) = v * i (r) (which, we note, is maintained throughout the encirclement process) plus the lack of overlap between the u's and v's on different vortices, guarantees that this is zero. So the total phase difference on encirclement is π and consequently that on exchange π /2, in agreement with Ivanov.
Note that for a non-encircling closed contour the evolution of delta at each of the vortices due to the other is zero, so we can repeat the argument and show that the phase difference due to such a process is zero. In the present (mean-field) approximation the π change on encirclement is a truly topological property.

The groundstate
We want to devote this first subsection to a question which, while it may or may not eventually turn out to be relevant to the main issue raised in this paper, needs to be resolved before we can claim a full understanding of a (p + ip) Fermi superfluid. The issue relates to the uniqueness of the many-body groundstate: In Sect. 2.1 we assumed that the groundstate wave function was given by eqn. (9), with the operatorĈ † given by (10); in words, we start from the vacuum and create on it N/2 Cooper pairs in a (p + ip) state. However, we note that the mean-field Hamiltonian (3) involves only two types of term, namely the particle number distribution n k and the "anomalous averages" a † k a † -k and its c.c. Thus, if we can find an alternative form of the GSWF which reproduces, exactly or approximately, the values of these two quantities, we will recover the exact or approximate value of the groundstate energy. Let us then proceed as follows (cf. ref. [11,Appendix 6A]): Denoting as usual the Fermi wave vector of the normal groundstate by k F the normal groundstate (Fermi sea) wave function by |FS and the phase of the parameter v k in eqn. (5) by -ϕ k , we first simply multiply the standard ansatz (9), by the overall factor k<k F exp{iϕ k } which is equivalent to the prescription where now we definễ and (remembering that we are working in the approximation of particle-hole symmetry, μ = E F ), In words: rather than starting from the vacuum and creating N/2 Cooper pairs, we instead start from the filled Fermi sea and "kick up" N + pairs out of the sea into the empty states above the Fermi surface. Since the ansatz (45) for the many-body wave function differs from (9) only by an overall phase factor, it represents the same state and must give the same values for any physical quantities, e.g. the (true) angular momentum. However; if we adopt it we can no longer use the mean-field Hamiltonian: since particle number is conserved separately for the C-and D-terms, such an approximation would neglect precisely those scattering processes which we are most interested in. This problem can be remedied by denoting the ansatz (45) as N (N + ) and writing the generalized ansatz where the coefficient q(N + ) is real and smoothly varying in magnitude of a range of N+ of the order of N( /E F ) (for present purposes we do not need to specify it more precisely than this). This not only allows the mean-field approximation, but gives a value for its groundstate expectation value which differs from that calculated for the ansatz (9) only by terms which vanish in the thermodynamic limit. This is equally true for any quantities which in the mean-field approximation are functions only of |v k | 2 and of the anomalous averages F k . However, it is essential to note that the ansatz (48) is NOT simply (9) multiplied by a phase factor; it is a physically different state, and thus may give different expectation values for physically meaningful quantities. Indeed, if we formally recalculate the "pseudoangular momentum" L ps defined by eqn. (8), we see that in the limit of particle-hole symmetry considered here it is exactly zero! (in the more general case it would be proportional to some power of N + -N -).
What of the behavior of the coordinate-space manybody wave function at large distances? At least prima facie, this should depend only on the small-k behavior of the quantities u k and v k . Unfortunately, if one tries to determine the answer from intuitive arguments, one reaches a dilemma: On the one hand, the difference between (48) and N (N + ) should be negligible in the limit in question, and since as we have seen the latter is just equivalent to (9), the result for the long-distance behavior of the ansatz (48) should be identical to that of (9), i.e. to the "Read-Green" ansatz. On the other hand, the k → 0 limit of (48), unlike that of (9), differs from that of the normal groundstate only by a factor which is negligible in the limit k → 0, and the long-distance behavior of the latter is quite different from that of (9)! At the time of writing we regard the resolution of this dilemma as an open question.
In conclusion, we note that the question discussed in this subsection is, at least prima facie, logically independent from that treated in the rest of Sect. 3; for purposes of discussing the latter, we shall assume the standard forms of the particle-conserving and -nonconserving many-body groundstates as used in Sect. 2.1.

Doubts about SBU(1)S
In particle physics it is usually held [17] that certain physical quantities are subject to a "superselection rule", that is, that coherent superpositions of states of a closed system corresponding to different values of one or more of these quantities are not physically meaningful. These quantities include, among others, lepton number and total fermion number mod 2. In the present context we will be interested in phenomena occurring in metals at room temperature or below, and moreover will assume for simplicity that all nuclei involved are immutable and bosonic: then both "lepton" and "fermion" in the above sentence can be replaced by "electron".
However, as we have seen, current approaches to the MF problem, and indeed to many other aspects of superconductivity, while they respect conservation of electron number mod 2, do not respect it "absolutely"; this is immediately clear from the form of eqn. (15), which is a superposition of many-body states differing by an even number of electrons. Is this physically meaningful? Let us first clarify a rather trivial point, which however seems to be capable of causing confusion: If we consider a closed box, then indeed the total number of electrons must, according to ref. [17], be a good quantum number. However, if as is often the case in a real-life experiment, the box has leads attached, then this condition clearly does not need to hold. Does this mean that the superselection rule can be violated? Certainly not: if we consider the reduced density matrix of the box alone (not including the leads) and write it in a basis where one of the indices is total particle number N, then in this situation, while different diagonal elements can be simultaneously nonzero (even for states of different number parity), if we believe ref. [18] as regards the combined box-plus leads system, then all matrix elements which are off-diagonal with respect to N must be identically equal to zero, because different N implies a different number of electrons in the leads, and tracing over the state of the latter then leads to that result. This principle is violated by the many-body state described by eqn. (15) and by any mixture of such states.
Despite the considerations of the last paragraph, there have been many attempts in the literature to justify the violation of the electron number superselection rule, i.e. the idea that the use of states of the form (15) is not just a mathematically convenient trick but has a real physical meaning. These attempts actually go back, in a slightly different but closely related context, to the more than 50-year-old work of Aharonov and Susskind [18], and have been advanced particularly energetically by the late P.W. Anderson (see e.g. ref. [19]); a favorite line of argument rests on the observation that the principle of conservation of electron number has, at least at first sight, a status no different from that of other superselection rules such as that of total angular momentum, and that in the latter case we are usually quite happy to violate it when discussing e.g. ferromagnetism. A thorough discussion of the validity or not of this analogy would eventually lead us back to some subtle issues in the quantum theory of measurement; fortunately, we do not need to go there in the present context, since our aim in this paper is not to convince the reader that it is necessary to build a rigorously number-conserving theory of superconductivity but rather to ask whether; if that should turn out to be indeed the case, it would affect the currently accepted lore concerning Majorana fermions. So from now on we will assume the necessity of number conservation.
One thing should be made clear from the start: To the extent that we are interested only in the properties of the condensate in the absence of any fermionic quasiparticles (i.e. in effect, only of even-number-parity manybody states), then it seems unlikely that this will make any substantial difference. The reason is that starting from a number-nonconserving state such as that described by eqn. (15), we can always perform the "Anderson trick" to produce from it a state with definite particle number. This trick works not just for the s-wave case for which it was originally written down, but also for e.g. the (p + ip) one: Write v k = |v k | exp i(ϕ k +χ ), denote the corresponding many-body state by (χ) and write which may be verified indeed to have definite (even) particle number N. The same trick works for odd-numberparity states: starting from a standard textbook Bogoliubov-quasiparticle state (with odd number parity), we again add a variable χ to the Cooper-pair part of it, integrate over χ and thereby produce a state with definite odd particle number and the same Bogoliubov particle. So in these cases it certainly seems reasonable to regard the supposed "spontaneous breaking" of the U(1) symmetry as merely a technical trick to expedite the calculations; this can be confirmed by noting that the original BCS procedure for solving the mean-field Hamiltonian (3) can be carried out with only moderately more labor in a completely numberconserving scheme (see [11, ch. 5]). The real difficulty arises only when one needs to relate states of even and odd number parity, i.e. to examine the process of creating Bogoliubov quasiparticles out of the condensate, and it is to this question that we turn in the next subsection.

Bogoliubov quasiparticles in the number conserving scheme
Consider first the case of a simple spatially uniform system, so that the total momentum is a good quantum number. Suppose that we insist that both the groundstate and all excited states of the superconductor which we wish to consider should be eigenstates of particle number. Then the obvious form of the even-number-parity groundstate for the system at rest is eqn. (9): including the normalization, with the constraint k |c k | 2 1+|c k | 2 = N/2 (51) and the normalization factor n given by Then what operator, acting on the even-parity groundstate, will create a Bogoliubov quasiparticle in, say, a state with momentum q? The textbook recipe (19), namely of course does so, but the resulting state is a superposition of states with different eigenvalues of particle number (N + 1 for the first term, N -1 for the second). Moreover, since these two states are mutually orthogonal, the Hermitian conjugate operator α q produces a state vector of nonzero length, i.e. it is not a pure annihilator.
It is tempting to dismiss these difficulties as a pseudoproblem, on the grounds that in real life we never operate on the (closed) experimental system with single fermion annihilation or creation operators, but only with bilinear combinations of them. So let's explore this point by imagining that we indeed apply to the number-conserving groundstate (50) a perturbation of the form (V qq a † q a q + H.c.), where we imagine, for notational simplicity only, that the coefficient V qq is real and symmetric.
In Appendix A we carry out such a calculation, with the conclusion that while either a consistently particlenonconserving (mean-field) calculation or a consistently particle-conserving one gives reliable answers, a "hybrid" calculation, in which the many-body states conserve particle number but the Bogoliubov operators do not, gives results which not only differ from the textbook ones but disagree with experiment. We therefore reach the conclusion that if we are going to take particle number conservation seriously for the even-parity groundstate then we must do so also for the odd-parity excited states (Bogoliubov quasiparticles). The obvious default way of doing this [9] is simply to modify the expression for α + q (eqn. (53) by leaving the u q term alone but associating with the υ q term an extra operator which creates two extra electrons with total momentum zero: where the operatorĈ † q , which for the moment we do not necessarily identify with theĈ † of eqn. (50), is of the generalized form where, again c q (k) is not necessarily to be identified with c k in eqn. (50), even apart from its normalization. We then see that if and only if c q (k) ≡ c k for k = q, -q the above difficulties are resolved, e.g. that the revised destruction operator α q , when applied to the groundstate, indeed creates a state vector of length zero (strictly speaking, of length N -1 , but we assume that this can be replaced by 0 in the thermodynamic limit). It should be emphasized that this maneuver, with the default choice for the operatorĈ † q ofĈ † except that the sum over k omits q and -q, is by no means novel in the theory of superconductivity. Indeed, some early texts such as [20] emphasize the necessity of making it, as do some important papers such as ref. [21]. However, over the ensuing 40odd years the point seems to have been gradually forgotten. We suspect that a major reason for this is that at least at first sight, when one is dealing with a simple BCS s-wave superconductor and a spatially uniform situation (which is what most textbooks understandably start with), the extra Cooper pair basically makes no difference-intuitively, it goes away "to infinity" and affects physical quantities only to order N -1 . Thus taking number conservation into account seems to modify the standard BCS results only trivially, and most recent papers just prefer to forget about it.
When one tries to generalize these considerations to a spatially inhomogeneous situation which needs to be handled by the BdG equations, things become a little more complicated. The analog of the procedure used above is to leave the terms in u(r) in the creation operators γ † i alone and to add to the v(r)-terms an extra operatorĈ † i , i.e. the definition of γ † i is now In this case it is somewhat less obvious that the correct choice ofĈ † i is the original even-parity Cooper pair oper-atorĈ † , that iŝ where K(r, r ) is the quantity defined in eqn. (12), but for the moment let us assume for simplicity that this is so. Then for the ordinary (Bogoliubov) quasiparticles the extra Cooper pair "runs off to infinity" just as in the s-wave case, and its addition would appear to have no real consequences.
Turning to the case of a Majorana fermion, which is the special case u 0 , v 0 of eqn. (56) with u 0 (r) = v * 0 (r), what are the consequences of replacing (34) by (56) (with as dis-cussedĈ † 0 =Ĉ † )? (1) The form of the BdG equations (35) for the coefficients u 0 (r) and v 0 (r) is not affected. (2) While the Majorana particle is no longer its own antiparticle (thus destroying the analogy with particle physics), it is still a quantum superposition of a zero-energy Bogoliubov particle and a pure annihilator, and it is still possible to combine two M.F.s to form a zero-energy electron just as in the mean-field approach.
(3) Since the extra Cooper pair "runs off to infinity" just as for the nonzero-energy Bogoliubov quasiparticles, the local undetectability of the M.F. is unaffected. (4) Braiding: see next section. What ifĈ † 0 is not equal toĈ † ? Then, (1) is still true (but one needs an extra equation to determine the form of K(r, r ). (2) However, the excitation which they describe is no longer simply a quantum superposition of a zero-energy Bogoliubov fermion and a pure annihilator, and it is in general no longer possible to combine two Majoranas to form a real Bogoliubov fermion. (3) Depending on the form and support of the added Cooper pair, the local undetectability of the M.F. may or not be destroyed. (generally speaking, it will be destroyed if the support is on the order of the pair radius or less). So, crudely speaking, it looks so far as if the quantum computing program will be unaffected by the corrections to mean-field theory if and only if the added Cooper pair is "the same" as the pre-existing ones. In the next section we examine whether or not this is true also as regards the braiding properties of the M.F.

Physical braiding properties of Majorana modes in the number-conserving approach
In this section we face head-on the question which has motivated this whole paper: does correct attention to the necessity of conserving total particle number invalidate some or all of the conclusions concerning the properties of Majorana fermions (or more accurately of the zeroenergy Dirac fermions constructed from them) which have been drawn using the standard BdG approach, in particular the results concerning their braiding properties? We should state at the outset that at the time of writing we have been unable to give a completely definitive answer to this question, however we will try to make some suggestive remarks. First, let us elaborate a little on the description, which was already given briefly in Sect. 2.5, of the specific experimental setup we are going to consider in this subsection. We consider a superconductor with (p + ip) pairing filling a 2D cylindrical (disk-shaped) container, and suppose that we are interested in the subspace of the Hilbert space in which it contains two vortices whose helicity is oppositely directed to that of the internal Cooper-pair state itself. We assume that the effective penetration depth λ eff is very large compared to the pair radius ξ (note that this condition is automatically fulfilled for a thin enough disk), and suppose that the core of one of the vortices is positioned at the origin R 1 = 0, while the other is carried around it at a core distance |R 2 | = R 2 which is ξ but both λ eff and the disk radius R d . Note that in the present context the ratio of λ eff to R d is mostly irrelevant; if it ever matters, we will choose it for convenience to be large, so that the superconductor is essentially equivalent to a neutral Fermi superfluid.
It will be necessary in the following to distinguish two different processes by which vortex 2 is carried around the origin and thus around vortex 1. (As we will see below, these two processes are actually best viewed as the two ends of a continuum.) In "process 1", which is the process implicitly considered in ref. [16] and hence Sect. 2.5, we keep the disk itself fixed but apply some source of external potential, for example an STM tip, which couples to the vortex core, and move it around, so that the Hamiltonian contains a term V ext (R 2 -R p (t)) where the potential V(X) has a "sufficiently deep" minimum at X = 0 and R p (t), the position of the probe, is given in plane polar coordinates by |R p (t)| = const. = R p , θ p (t) = ωt. In "process 2" the whole disk, including the confining walls (and if relevant the crystal lattice) is rotated at angular velocity ω. In either case the process is assumed to be adiabatic, meaning that if v denotes a characteristic velocity of the system such as the normal-state Fermi velocity or the speed of hydrodynamic sound, the angular velocity is small compared not just to v/R 2 but to v/R d . Under these conditions we assume that to a sufficient approximation the angular position θ 0 (t) of the core of vortex 2 is identical to θ p (t), i.e. is given by ωt.
The rotation starts at time 0 and ends at T = 2π/ω, corresponding to a single encirclement.
As stated in Sect. 2.5, the problem of interest (which as there we shall call the "Ivanov problem") is to determine the difference in the phases accumulated under the encirclement process by the even-and odd-parity groundstates (without and with the zero-energy Bogoliubov quasiparticle respectively). We will call this difference η.
Let us start with some statements which we believe to be unproblematic: (1) We first consider the even-parity groundstate in equilibrium (θ 0 fixed, say for convenience at 0) and write the anomalous average (order parameter, OP) F(r, r ) ≡ ψ (r )ψ(r) as a function F(R, ρ) of the COM variable R and the relative coordinate ρ; note that F must be antisymmetric in ρ but can have a single-valued but otherwise arbitrary dependence on R. We shall assume that everywhere in the region R 2 R λ eff the OP is approximately a product of the standard (p + ip) dependence on ρ with helicity directed along the negative z-axis, and a function of the COM coordinate R which is approximately constant as a function of |R| but rotates clockwise twice as a function of the angle theta with respect to the cylinder origin: where F 0 is a constant (possibly complex, see below) amplitude and f p+ip is the standard coordinate-space OP for a p + ip pairing state as discussed in Sect. 2.1. For R within a distance ≤ ξ of either the origin or R 2 , the OP is more complicated, but has the property that when R rotates around either the origin or R 2 its phase rotates through 2π . Note that under these conditions the total (COM plus internal) angular momentum is -N /2 to within corrections of relative order (R 2 /R d ) 2 (a qualification we shall omit for the following discussion). (However, we should bear in mind that this is only true to the extent that we assume the "standard" form of f p+ip , cf. Sect. 3.1.) We note that the numberconserving many-body wave function has the form, up to normalization, where the operatorĈ † is given by the expression The function K(r, r ) entering this expression is not identical to F(r, r ), but has the same topological dependence on R and ρ. (For a detailed discussion of the two quantities, see e.g. ref. [22], Appendix A: cf. also eqns. (22)-(27).) (2) Next, we discuss the addition of a single extra Cooper pair to the even-parity groundstate. This may be achieved by applying the operator where for the moment we do not assume thatK(r, r ) = K(r, r ). The extra angular momentum added is given by and some straightforward algebra yields as intuitively expected.
(3) Turning to the odd-parity state, we note that since that the BdG equations remain valid for the calculation of the coefficients u(r) and v(r), all we need to know to calculate either the values of these in equilibrium or the monodromy phase due to encirclement is the value (or change due to encirclement) of the OP (and hence of the gap) in the regions of the vortices, that is within ∼ ξ of the origin or of R 2 (we do not need, at this stage of the argument, to know the evolution of these quantities over the whole of the disk, or even over the whole region of radius R 2 around the origin). In particular, the BdG equations plus the condition of self-conjugacy (which must be preserved during the rotation) together imply that the monodromy phase is exactly half of the phase acquired by (or equivalently by F).
(4) If and only if the correct form of the wave function during encirclement is the "fully rotated" one of eqn. (65) below, then there is a simple relation [5] between the angular momentum of a given state and the Berry phase ϕ B due to encirclement, and thus also between their differences: All the above statements are true independently of the mean-field approximation, and we believe are not controversial.
Let us now turn to the question: How is the many-body wave function, and hence the OP changed, by processes 1 or 2 respectively? Here our arguments are plausible rather than rigorous. We actually start with process 2 which is conceptually rather simpler: (5) Process 2: Since in this ansatz "everything" is rotated, the obvious ansatz for (t) N and hence for F(R, ρ : t) is 4 where θ R , θ ρ are the polar angles respectively of R and ρ.
In words, both the COM and the relative coordinate are rotated through ωt.
With this ansatz, it is clear that under encirclement (of one vortex by the other) the OP at each vortex due to the other is rotated through 0 (+2π for R and -2π for ρ).
(6) Process 1: this is a little trickier, and indeed the correct answer may depend in detail on factors like the "strength" of the rugosities on the disk edges etc. However, in the limit that R 2 /R d 1 a plausible ansatz is that at least within a distance ∼ R 2 of the origin, or in words: the COM coordinate is rotated but the relative coordinate is not. If this is the case, then the OP and hence is rotated through -2π under encirclement. Armed with the results (1)-(6), we are now in a position to discuss the Ivanov problem, both within the mean-field scheme and with the constraint of particle number conservation.
First, we revisit the mean-field scheme of ref. [16] and Sect. 2.5, in which γ † i involves only fermion operators.
(a) Process 1: According to (6), at the position of the vortices the only change of the OP is a rotation of by 2π , so as detailed in Sect. 2.5, the monodromy phase is π and the Berry phase 0. Thus the total phase difference is π in agreement with Ivanov. (b) Process 2: Now, by (5), the rotation of delta at the positions of the vortices is 2π , but that of the internal OP is -2π , so the monodromy phase is zero (note not 4π !) and the Berry phase is also zero. Thus the total phase difference is zero, and that on exchange also zero, in contradiction to Ivanov. Now we turn to the crunch: the effect of demanding respect for particle number conservation. The only effect on the even-parity groundstate is that we must now use the particle-number conserving version N . What about the individual Majoranas? ln order to conserve particle number, we must now rewrite the Majorana creation operator γ † i , for the MF which resides on vortex i in the form where for the momentĈ † i is of the general form eqn. (61). We note that in the general case this maneuver, apart from destroying the self-conjugacy of the M.F. 's and hence the analogy with particle physics, has the two immediate consequences listed at the end of Sect. 3.3.
(2) It invalidates the relation between the M.F. and the zero-energy Bogoliubov quasiparticle. (3) It creates an extra particle density and thus destroys the local undetectability of the M.F. Since the "default" optionĈ † i =Ĉ † , whereĈ † is the operator (eqn. (60)) which enters the expression (eqn. (59)) for the even-parity groundstate, does not have the consequences (2) and (3), let us for the moment make that ansatz and return later to the more general case.
So, back to the Ivanov problem. We first note that since the gap is sensitive only to order N -1 to the presence or absence of the extra pair, in both processes the arguments concerning its evolution under encirclement are unchanged and since so are the BdG equations, the conclusion concerning the monodromy phase is identical to that in the mean-field case, namely π for process 1, 0 for process 2. However, the considerations concerning the Berry phase are changed: for process 2 the extra Cooper pair adds an extra angular momentum 1 to within terms of relative order (R 2 /R d ) 2 , but this is associated only with "half " of the Bogoliubov quasiparticle, so L is 1/2 and the Berry phase correspondingly π , giving an overall phase difference η = π in agreement with Ivanov.
For process l, the result depends on how the wave function is distorted in the bulk of the system by the rotation: if this is either negligibly small, or has the form (eqn. (67)), then the Berry phase is 0 to within terms of relative order (R 2 /R d ) 2 and hence we again recover agreement with Ivanov. If however the assumption in the last sentence fails (or if R 2 is not negligibly small compared to R d ), then to obtain the phase difference we really need a proper solution of the time-dependent Schrödinger equation for the problem of the moving probe, which at the time of writing we do not have. However, it seems at least a nonnegligible possibility that if the experimental encirclement process is indeed either process 1 or something intermediate between processes 1 and 2, the total encirclement phase may not be not a "nice" number (i.e. a multiple of π ) at all, but may be a complicated function of the system parameters. Needless to say, should this turn out to be the case it would be bad news for topological quantum computation using 2D (p + ip) Fermi superfluids.
The problem becomes an order of magnitude worse if we consider the possibility that the operatorĈ † i which adds the extra Cooper pair is not the "default" optionĈ † , but contains some antisymmetric but otherwise arbitrary function K(r, r ).
One prima facie attractive option is to restrict the support of K(r, r ) to a distance of order ξ from R i (i.e. to the vicinity of vortex i). Clearly, such an ansatz would give results much different from the ones obtained above; indeed, the whole TQC operation in this system would not get off the ground.
Thus, it becomes of the utmost importance to determine the nature of the extra Cooper pair added to the v-component of a Majorana fermion. Unfortunately, this looks like a complicated problem, since the moment we fail to assume the "default" form ofĈ † i we need to supplement the BdG equations for the u-and v-components with an extra equation, so that it is not even obvious that Majorana solutions will exist at all. And intuition does not really seem to be of much help: at first sight it seems that "common sense" would imply that a Cooper pair created to keep the auditing straight for a MF on vortex 1 should not "know" that it needs to have orbital angular momentum -1 about the distant vortex 2 (an argument which, if true, would favor the "localization-by-ξ " hypothesis); on the other hand, what distinguishes the so-called topological superfluids is precisely this kind of long-range "knowledge" (entanglement)! If the real problem looks intractable, one tactic may be to try to solve a simpler model problem which has some of the right features. In ref. [6] we attempt such a solution for the problem of a mesoscopic Josephson junction with an odd number of particles (hence total spin 1/2) in a localizing magnetic field, and find that at least a certain amount of localization (of order unity in the thermodynamic limit) of the extra Cooper pair accompanying the v-component of the fermion does take place (suggesting the "localizing" answer). However, the problem is sufficiently different from ours that it may not be legitimate to transpose this answer.

Conclusion
As emphasized in the Introduction, this paper has been devoted to a rather specific question: does the simplest modification of the "standard" theoretical approach to Majorana fermions in a 2D p + ip Fermi superfluid which guarantees conservation of overall particle number modify their properties, and in particular their behavior under "physical" braiding, in any essential way? We have seen that the answer to this question may depend to some extent on the precise definition of "physical braiding", and in particular whether it is instantiated by process 1, process 2 or something intermediate between these two processes. In the case of process 2, we have seen that a very plausible answer is that while the mean-field calculation gives an answer different from the mean-field ("Ivanov") one, the simplest number-conserving modification restores the latter. In the case of process 1 (the process implicitly considered in ref. [16], and probably most likely to be used in a real experiment) it appears probable that the result of the braiding process will depend on the details of the experimental geometry, with the mean-field result emerging only in a particular limit.
However, one general point needs to be noted: although we have spent a lot of time, in Sects. 2.5 and 4, discussing the difference in the phase acquired by the even-and oddparity states under the encirclement operation, it is not clear that this is itself a measurable quantity; in any real-life attempt at topological quantum computing, we are likely to be braiding, i.e. operating within a subspace of constant total particle number N, and the solution of the "Ivanov problem" is at best an intermediate step in the calculation of such a process. It seems entirely possible that while the formal arguments of (e.g.) ref. [16] need to be substantially modified to accommodate number conservation, at the end of the day the TQC protocol still works as advertised. As far as we know, the existing literature has not addressed this question (cf. Appendix B), and the constraints of time and space require that we leave this important question for another day.
While these (tentative) conclusions may be of some interest in the context of the specific question studied, from the overall standpoint of the topological quantum computing enterprise this question may well by now look rather parochial. It seems that future progress in that field, even in so far as it relies on p + ip (-like) Fermi superfluids at all, may well prefer the 1D version (and perhaps the use of measurement-based braiding), and the differences in the physical (as distinct from the algorithmic) aspects of the problem both between 1D and 2D systems, and between the physical and measurement-based implementations of the braiding procedure, are sufficient that in our view it would be risky at present to draw conclusions too firmly in either direction. Whether this situation will remain true as the subject develops remains to be seen.

Appendix A
Effect of particle nonconservation in intermediate states on particle conserving processes.
In this appendix we consider the results of treating a single-particle perturbation of the form ( qq V qq a † q a q + H.c.) acting on a Cooper-paired system (a) if both the groundstate and the excited states are given the "textbook" (particle-nonconserving) treatment (b) if the groundstate is required to be an eigenstate of particle number but the Bogoliubov quasiparticles are given the "textbook" treatment, and (c) if particle conservation is strictly required for both the groundstate and the Bogoliubov excitations. For simplicity we take V qq to be simply a constant V 0 , so that time-reversed terms enter with the same sign. We will consider nonzero as well as zero temperature, so that the density matrix may include "unpaired" states, and focus on the single term of the form α + q α q which describes the scattering of a quasiparticle from state q to state q.
The answer to (a) follows immediately from textbook material: since in the BCS approach the density matrix is a product of sub-density matrices referring to different pair states (k, -k), the calculation can be confined to the subspaces associated with (q, -q) and (q , -q ), the Bogoliubov transformation takes its standard form and we simply get for this term the expression (up to complex conjugations) V 0 u q u * q -v q v * q α † q α q .
Next consider question (b). ln writing the many-body state we now need, in order to guarantee overall particle number conservation to take into account explicitly the states of all the pairs other than q, q . Denote the number of pairs in such other states by M qq (it is actually a complicated state vector but this does not concern us at present). Then, for example, the initial energy eigenstate corresponding to a ground pair in state q and a broken pair in state q should be written not (as in the BCS approach) |10 q |GS q times the groundstate of the "other" pairs, but rather in = |1, 0 q · u q |00 q |M qq = N -1 + v q |11 q |M qq = N -3 .
Now if we first blindly apply the textbook Bogoliubov transformation (19) to the perturbation and then apply it to this state, we find that the resulting coefficient of the final state is not as above (u q u * q -v q v * q ) but rather u q u q |M qq = N -1 -v q v * q |M qq = N -3 .
Consequently, since states with different M qq , are mutually orthogonal, the squared matrix elements between the initial and final states which in the standard notation would be α † q and α q times "something irrelevant" will contain not |(u * q u q -v q v * q )| 2 but rather u q u * q |M qq = N -1 -v q v * q |M qq = N -3 just the diagonal terms |u q u * q | 2 +|v q v * q | 2 . In general this will lead to different expressions for the effect of various probes in the superconducting state, and thus (since the standard expressions are well verified experimentally) we must conclude that ansatz (b) is excluded.
Finally, we see that as noted in the text, it is very straightforward to fix up ansatz (b) so as to avoid the above contradiction with experiment: We simply associate with each factor v q in the Bogoliubov transformation an extra operator C † which creates an extra Cooper pair; where for the moment C † is just the quantity which enters the groundstate wave function in the particle-conserving approach. It is easily seen that this gives back the standard result, coherence factors and all.
The above material may no doubt seem rather trivial, but the moral is that if we are going to insist on enforcing particle number conservation for the even-and odd-parity states separately then we must also do so for their relation, and that failing to do so may lead to physically incorrect results even for processes which, taken overall, conserve particle number. This remark becomes less trivial when applied to calculations involving Majorana fermions

Appendix B: Comments on some recent papers
While as noted in the Introduction the discussion of the role of overall particle number conservation in the theory of superconductivity is as old as the BCS theory itself, its examination in the context of topological superconductors is relatively recent. We may classify existing papers into three categories: (1) Papers which consider operations such as (physical) braiding in 2D p + ip Fermi superfluids and pay some attention to number conservation. (2) Papers which address the parity question for a number-conserving 2D p + ip topological superfluid. (3) Papers which address the relation between different-parity states (and in some cases more) in a number-conserving 1D model. We discuss representatives of each of these classes in order: (1) The paper closest to our own work in this category of which we are aware is ref. [23], which while it shows that a particular (rather exotic) number-conserving Hamiltonian can generate exact mean-field groundstates, then carries out the braiding operation within the even-parity sector but not for definite N. We are unaware of any paper in the existing literature which does the latter (for 2D). (2) The most ambitious member of this class to date seems to be ref. [24]; this deals with both 1D and 2D systems, and at this point we discuss only the results