Probabilistic Q-function distributions in fermionic phase-space

We obtain a positive probability distribution or Q-function for an arbitrary fermionic many-body system. This is different to previous Q-function proposals, which were either restricted to a subspace of the overall Hilbert space, or used Grassmann methods that do not give probabilities. The fermionic Q-function obtained here is constructed using normally ordered Gaussian operators, which include both non-interacting thermal density matrices and BCS states. We prove that the Q-function exists for any density matrix, is real and positive, and has moments that correspond to Fermi operator moments. It is defined on a finite symmetric phase-space equivalent to the space of real, antisymmetric matrices. This has the natural SO(2M) symmetry expected for Majorana fermion operators. We show that there is a natural physical interpretation of the Q-function: it is the relative probability for observing a given Gaussian density matrix. The distribution has a uniform probability across the space at infinite temperature, while for pure states it has a maximum value on the phase-space boundary. The advantage of probabilistic representations is that they can be used for computational sampling without a sign problem.


I. INTRODUCTION
Phase-space representations of the quantum mechanical density operator were introduced by Wigner [1] and Moyal [2]. Such mappings from quantum mechanics to phase-space have many physical applications, especially in coherence theory [3]. Positive phase-space distributions or Q-functions, introduced by Husimi [4], give mappings which allow quantum mechanical observables to be calculated using probabilistic sampling.
This paper introduces a unique, positive fermionic phase-space representation. This has the same useful properties as other Q-function methods which have been used for bosonic and spin Hilbert spaces. The Q-function derived here is a unique positive distribution, defined for all density matrices, and applicable to fermionic manybody systems. Since it is probabilistic, there are no sign problems when sampling the distribution.
Examples of this general approach include the bosonic Q-function [4], the SU(2) Q-function [5,6], the positive P-function [7], and the general positive distributions over Gaussian operators [8][9][10]. It is known that the sampling properties of Q-function methods scale favorably in the limit of large systems, which allows high order correlation functions to be computed for large spin systems [11]. Q-functions using atomic coherent states have also been used to study the dynamics of superfluorescence [12][13][14]. Another application of this type of distribution, based on coherent states for spin and oscillator states, has been used to study quantum dynamics of the Dicke superfluorescence model [15] and thermalization processes [16].
We introduce a general method for constructing such probability distributions, which is not restricted to Fermi systems. Our definition is based on the expectation value of a hermitian operator basis. In the present case, this is the basis of Gaussian fermionic operators [9,10,17]. Importantly, we prove that these provide a resolution of unity [18], and have simple differential identities for operator moments. This is fundamental to constructing a useful probabilistic representation.
A Q-function provides a useful method for visualizing and understanding coherence and correlations. It also has potential applications as a computational tool. Since the distribution is always probabilistic, there are no intrinsic sign problems with sampling or Monte Carlo methods. We will treat detailed applications elsewhere.
To achieve this result, we require the representation to satisfy the following three properties, all found with the Husimi function: 1. It exists uniquely for any quantum density-matrix.

2.
It is a positive probability distribution.

Observables are moments of the distribution.
Can we satisfy these conditions in the case of fermions?
Here, we obtain a fermionic Q-function which satisfies the three requirements given above, and show how this fits into a general picture which includes the well-known Husimi function. Our method makes use of a complete set of positive-definite normally ordered operators as a basis [10,17,18], and extends the conceptual basis of Q-function methods. The phase-space is a bounded domain of real antisymmetric matrices. We show how the theory of matrix polar coordinates and Riemannian measures on symmetric spaces provides a natural mathematical framework for this probability. Physically, we make use of the Class D symmetries introduced by Altland and Zirnbauer [19] to treat normal-superconducting interfaces. We focus on this case, as it is the most general, while bearing in mind that similar results will follow in each of the four symmetry classes of this type. This paper is organized as follows. In the next section we introduce our general approach of defining a Qfunction. In section (III) we describe the general normally ordered Gaussian operators and their properties. In section (IV) a completeness proof is given, based on matrix polar coordinate theory. These results are used to define the Q-function in section (V), together with identities for observables and moments. In section (VI), we give explicit expressions of the formalism developed here for Gaussian density operators. Section (VII) gives a summary of our results and conclusions. Finally, properties of the unitary transformations that we use are proved in Appendix (A) and the inner product of the general Gaussian operators is calculated in Appendix (B).

II. GENERALIZED Q-FUNCTION PROPERTIES
We first introduce a general definition of a positive Q-function representation for any Hilbert space. The approach given here is not just applicable to fermions. It is therefore useful to understand the abstract concepts first. Our approach differs in an essential way from that of Husimi [4] and others [5,20]. Rather than constructing the phase-space representation using coherent states, we employ a hermitian operator basis. This operator basis can be thought as a basis of coherent operators, rather than coherent states. It is closely related to the more general idea of a Stratonovich-Weyl mapping [21,22].

A. Q-function definition
Suppose we have a positive definite, hermitian operator basisΛ λ defined in a Hilbert space H of quantum mechanical operators, where λ is a vector in the phasespace domain D. We require the following completeness property: the identity operatorÎ of the Hilbert space can be resolved as an integral over the phase-space, so that DΛ λ dµ λ =Î .
(2.1) This is called a resolution of unity. Here dµ λ is an associated integration measure on the phase-space. A generalized Q-function is defined as the inner product of the density matrixρ with the operator basis: With this definition, condition (1) is automatically satisfied. We will show below that condition (2) is also satisfied. Hence, a Q-function exists for any density matrix in a Hilbert space which has a continuous, positive resolution of the identity operator, proved suitable identities are found that satisfy condition (3).
There is a direct physical interpretation. The opera-torsΛ λ are observables, and the distribution simply gives their measured expectation values. In the next sections, we demonstrate how the conditions given above can be satisfied in the case of fermions.

B. Probability distribution
To demonstrate that the generalized Q-function defined above has the properties of a probability distribution, as required for condition (2), we must prove that it is positive and normalized for any density matrix.
a. Positivity: Since the density matrixρ is hermitian and positive definite, it must have a diagonal Schmidt decomposition in an orthogonal basis |n , as: ρ = n ρ n |n n|, where ρ n ≥ 0. SinceΛ λ is positive definite by assumption, n|Λ (λ) |n ≥ 0, so that Thus, the resulting distribution function is positivesemidefinite, as required.
b. Normalization: From the definition given in Eq. (2.1), it is clear that the trace of any operator can be expressed as an integral over the phase-space, since: Tr Ô =ˆTr ÔΛ λ dµ λ . (2.4) ChoosingÔ =ρ, it follows that the distribution is normalized to unity, since from Eq (2.2):

C. Observables and moments
Next, we wish to satisfy condition (3) for the general Q-function case. This requires the evaluation of observables in the form of Tr ρÔ n = Ô n , and it leads to nontrivial requirements on the basis set. Since the eigenvalue methods used for the Husimi Q-function are not always available, we look for a more general approach. Proving that these requirements are satisfied in the case of fermions will require an understanding of the differential properties of the fermionic Gaussian operators.
We suppose that there are a complete set of identities that allow all operator moments of physical interestÔ n to be mapped into differential operators, so that: Here the second equation follows from the first, together with the assumption that the basis set is hermitian. We use the convention that in D n ∂ λ , λ , all differential operators are ordered to the left of any functions of λ. It then follows using the resolution of the identity, Eq (2.1), that any observable in the form of an operator moment can be represented as: Finally, provided D n ∂ λ , λ can be transformed via integration of the differentials into O n λ = D n 0, λ , with vanishing boundary terms, it follows that: If the differential mappings exist, and the distribution Q λ allows partial integration, observables can be calculated as a probabilistic distribution of moments over λ. It is important to choose a phase-space mapping such that O n λ is efficiently computable.
The set of observablesÔ n needs to include only the physically relevant moments. If the Hamiltonian has conservation laws, we are usually not interested in their invariant dynamics. This allows one to reduce both the Hilbert space and phase-space dimensionality. Dimension reduction proves useful in the case of the ground state properties of the fermionic Hubbard model, which was numerically solved using fermionic P-function methods using both translational and number-conservation symmetries [23].

D. Bosonic Q-function
Before examining the fermionic case, we now show that this definition includes the well-known bosonic Qfunction [4]. In this case the relevant operator basis is a normalized set of bosonic coherent state projectors | α α|, which are also used in the Glauber-Sudarshan P-representation [24,25]. These are defined in terms of the normalized bosonic coherent states | α , which are eigenstates of the bosonic annihilation operators, a = (â 1 , . . .â M ) with eigenvalues α = (α 1 , . . . α M ): These match our definition: they give a resolution of the identity, since they have the property [24] that 1 π Mˆ| α α| d 2M α =Î . (2.10) The relevant phase-space of λ is the M −dimensional complex space C M , and the measure is the standard Euclidean volume measure of d 2M α = j dα (y) . From the definition in Eq (2.2) above, one can express the Husimi Q-function as a unique positive distribution: In this case, the mappings to observables are also straightforward. For antinormally ordered operators of formÔ mn = â mi i â †ni i , the corresponding function on phase-space is a c-number moment O mn ( α) = α mi i α * ni i , so that This result is obtained, as usual, from the application of the eigenvalue equation for the annihilation operators and the definition of the Q-function.

III. FERMIONIC GAUSSIAN OPERATORS
There have been several previous approaches that have reproduced some, but not all of the properties described in the Introduction. One approach [26] was to directly introduce an SU (2) based fermion coherent state. This satisfied (1) and (2), but not (3). Subsequently, a Qfunction based on fermionic coherent state projectors was introduced using U (N ) Lie group methods [6,[27][28][29][30]. This satisfied (2), but not (1) and (3), since the projectors are not a complete basis, and observables were not derived. A third approach due to Cahill and Glauber [31], used Grassmann coherent states. This approach satisfies (1) and (3) but not (2), as a Q-function defined this way is Grassmann valued, and therefore neither real nor positive.
Here we define the fermionic Q-function using the method given above, with fermionic Gaussian operators [17]as a basis. These are also utilized in the complementary fermionic P-function representations [10], which are non-unique mappings to a phase-space of larger dimension. These have been utilized for evaluating the ground state of the fermionic Hubbard model, via Monte-Carlo techniques [23]. Other applications of these fermionic P-functions are to the quantum dynamics of Fermi systems, like molecular dissociation [32,33], and the linear entropy in a quantum phase space [34].

A. Quadratic Hamiltonians
Consider a fermionic system composed of M modes. We defineâ as a vector of M annihilation operators and a † as vector of M creation operators, whereâ i andâ † j obey the fermionic anticommutation relations: An extended vector of all 2M operators is written aŝ a = â T ,â † T , while the corresponding adjoint vector iŝ a † = â † ,â T = â † 1 , . . . ,â † M ,â 1 , . . . ,â M . We denote 2M × 2M matrices as H, and M × M matrices as h.
The most general quadratic form of the fermion fields, expanded in mode operators, is just the well-known Bogoliubov-de Gennes Hamiltonian: This can written in a compact form asĤ = 1 2â † Hâ, where the 2M × 2M matrix H is defined as: In order forĤ to be hermitian, one has the fundamental requirement that h = h † and that ∆ = −∆ T . These conditions can be written [19] as conditions on H, where: A Hamiltonian with ∆ = 0 is number-conserving, and describes a non-interacting Fermi gas. With ∆ = 0, this becomes a phase-dependent form appropriate for describing normal-superconductor interfaces, where Cooper pairs can tunnel through a barrier, and is also useful in mean-field treatments of superconductors.

B. Gaussian operators
Gaussian operators arise in many physical contexts. In analogy with Gaussian distributions in probability theory, these are defined to be exponentials of quadratics in the field operators [9]. Here, we will introduce a unit trace, hermitian Gaussian operator proportional to exp −βĤ . Although these also play the role of freefield canonical density matrices, the states that can be represented are completely general.
Our motivation is similar to Glauber's [24] use of bosonic coherent state projectors, which are also Gaussian operators [8]. Just as with bosons, fermion states with phase-dependent terms do not physically occur in isolated systems. However, they describe coherence properties in the simplest possible way, and provide a complete basis which can be used to represent any density matrix.
It is convenient to parametrize such Gaussian operators as unit-trace normally ordered forms, defined as [10,17] whereΛ u is a normally-ordered but un-normalized Gaussian operator,Λ andĪ is a diagonal matrix given by: Here 0 and I are the M × M zero and identity matrices, respectively. Normal ordering is denoted by : . . . :, hence is the covariance matrix, which has an identical symmetry to H. In terms of an M × M hermitian matrix n and a complex antisymmetric matrix m, one can write: In physical terms, the Gaussian operatorΛ is simply the normally ordered density operator of a finite temperature, noninteracting Fermi gas with Hamiltonian H, where n is the normal Green's function and m the anomalous Greens' function. These density operators arise in the theory of non-interacting Fermi gases and BCS superconductors. The relationship between the σ and H matrices is known [35], as is the fact that the symmetry properties of σ and H are the same. However, for our purposes, only the normally ordered covariance matrix σ is important, as it will define the Q-function phase-space.

C. Differential identities
The normally ordered Gaussian operators have the advantage of having simple differential identities, which will allow us to obtain the Q-function observables. These correspond to the action of the extended creation and annihilation operators on the Gaussian basis [10,17]: : Hereσ ≡Ī − σ. We use nested ordering {: . . . :} with the convention that external orderings do not change orderings inside internal brackets, and the ordering ofΛis invariant. When calculating matrix derivatives, we use the convention that: ∂/∂σ αβ ≡ ∂/∂σ βα [17]. Due to matrix symmetry, one finds that

D. Group properties and transformations
In order to prove our results in later sections, it is useful to relate the Gaussian operators to fundamental results on the symmetry properties of physical systems, and the concept of a symmetric space.
Dyson's 'threefold way' [36] classified the symmetry groups of physical systems into real, complex and quaternion types, corresponding to the Weyl classical groups [22]: orthogonal, unitary and symplectic. The physical meaning of these three classes relates to behavior under time reversal and spin rotation. Additional, nonstandard symmetry groups have been found since, making ten in all. The work of Altland and Zirnbauer [19] considered symmetry classes arising in the case of quadratic Hamiltonians for coupled superconductornormal fermionic systems. These are just the transformations of the Gaussian fermionic operators defined above.
The most general case has neither time-reversal nor spin-rotation invariance. This case, known as class D symmetry, allows one to treat an arbitrary Fermi system.
Many physical systems have a higher degree of symmetry than this. This can be used to reduce the phase-space dimension, which has practical advantages. For example, number conservation and translational symmetry are utilized in Imada's analysis of the Hubbard model ground state [23], which uses Gaussian operator methods. Such dimension reduction methods can used for Q-functions also, but we do not treat them here for simplicity.
All Class D matrices have the symmetry indicated by Eq (3.4). To obtain expressions which treat particles and holes on an equal basis, we introduce the zeta matrices: which we will call 'stretched' variables. These have an identical class D symmetry.
In order to understand their relationship with the classical symmetry groups, we can transform the Fermi operators to an hermitian Majorana fermion basis,ẑ = U 0â . This allows us to introduce a corresponding matrix mapping that preserves quadratic forms: where the unitary matrix U 0 is defined as: (3.14) The Class D symmetry properties then become X = X * = −X T , so these are 2M × 2M real antisymmetric matrices. Like the matrices iζ, they form an so(2M ) Lie algebra. This implies that the Gaussian operatorsΛ are a representation of the SO(2M ) Lie group [37]. While we focus mainly on the usual Fermi operators here, clearly there is an equivalent approach using hermitian Majorana operators and the real antisymmetric X matrices.
An important consequence of this group theoretic correspondence [19] is that there is a unitary matrix U that diagonalizes any hermitian matrix ζ, while retaining the symmetry of Eq (3.4). This unitary matrix is also a member of the SO(2M ) Lie group, with the property that This requirement means that the unitary transformation b = Uâ preserves the Fermi commutators [38].
The preservation of Fermi commutators on the extended vector of operators is essential to the unitary transformations in this paper. It is used in Appendix A to show that a normally-ordered Gaussian operator remains normally-ordered in the new variables, after a unitary transformation of both the operators and the covariance matrix.

E. Positivity
WhileΛ is hermitian, is it also positive-definite? Because not all Gaussian operatorsΛ of the form given in (3.7) are positive definite, only a finite domain in the phase-space of ζ variables leads to a physical density matrix with positive eigenvalues. To find a necessary and sufficient restriction, we use the fact that the corresponding covariance σ matrix is hermitian, and so is diagonalizable using the unitary transformation of Eq (3.15) to new operatorsb = Uâ, leaving the commutation relations invariant.
We show in the Appendix that this diagonalizing transformation can be applied inside the normally ordered symbols, and leaves invariant the Class D reflection symmetry on the diagonals. After transforming ζ, we obtain where n ′ is diagonal. On re-ordering the anti-normal terms in the exponential, which requires a sign-change owing to the definition of fermionic normal ordering, the un-normalized Gaussian becomes: (3.17) Noting that the possible eigenvalues ofb † jb j are 0, 1, we see that the condition (1 − n ′ ) n ′ > 0 is both necessary and sufficient for positivity ofΛ u and henceΛ. Since a unitary transformation does not change the positivity of a matrix, this can be re-written as an equivalent condition on ζ, defining the phase-space domain D as: The stretched ζ matrices therefore have an eigenvalue range from (−1, 1). This domain corresponds to the positive definite hermitian operators of interest.

F. Classical domains and symmetric spaces
Since the class D symmetric covariance matrices can be transformed into real antisymmetric 2M × 2M matrices using Eq (3.13), the hermitian Gaussian operators can also be regarded as having a real phase space of M (2M − 1) dimensions. The boundary of this space, equivalent to the domain D of Eq (3.18), is defined by the requirement: This is the real subspace of the irreducible homogeneous bounded symmetric domain R III of complex skewsymmetric matrices [39], called a classical domain in the theory of matrix polar coordinates [40]. More generally, all spaces of this type, with Riemannian measure, have a one-to-one relationship with the ten simple Lie groups that describe physical symmetries. These are called the symmetric spaces [41].
Geometrically, we note the following property in the present case: Thus, every individual element of X is bounded, since |X ij | < 1, and every row and column is bounded by a corresponding hyperspherical shell. This shows that the Class D group symmetry properties have a natural correspondence to a phase-space with a finite boundary. We note here that there are four similar types of symmetric space of this general class, with differing symmetry properties [19]. In principle, any of these can be used to coonstruct a different type of fermionic Q-function that is appropriate to the relevant symmetry group.

IV. RESOLUTION OF UNITY
In order to obtain the Q-function representation, the next step is to obtain a resolution of unity for the normally-ordered fermionic Gaussian operators expressed in terms of the ζ matrices. We follow a similar procedure to the parity argument used previously to obtain the resolution of unity for the fermionic Gaussian operators [18]. The main difference is that the present identity requires an integral over a finite domain, due to our choice of a normally ordered basis set.

A. Riemannian volume
To carry out integration on the phase-space, it is necessary to obtain a volume measure. Here we follow the original approach of Hurwitz [42,43] and Hua [40], by using measures defined relative to an invariant metric This gives a unified invariant measure over both the real and complex matrices described above. Given a distance metric, ds 2 = g ij dx i dx j , the corresponding Riemannian volume measure is dµ = |g| dx i [44].
To start with, we define: which are the Euclidean measures for the independent elements of the real anti-symmetric matrix X, and the class D complex matrices respectively. The volume elements can now be computed. In the antisymmetric case, the metric and corresponding Riemannian measure are For Class D hermitian matrices the metric and Riemannian measure are larger, since: Thus, we have now defined the phase-space and a volume measure for the Q-function operator basis.

B. Matrix polar coordinates
Integration over the phase-space of matrix variables is simplified by using matrix polar coordinates, which are commonly used in random matrix theory [19,40,45]. Since iζ belongs to a Lie algebra, and ζ is hermitian, it can be diagonalized, see Section (III D): Here U is an element of the SO (2M ) group and ζ D = diag (ζ, −ζ). The eigenvalues must be in the range −1 < ζ j < 1, since the domain is such that 1 − ζ 2 > 0. The Jacobian for the transformation from the coordinate ζ to polar coordinates ζ, U , is given by [19,45]: Here we have used the matrix polar coordinate [19] measure for class D symmetry, where dζ = M j=1 dζ j and ∆ ζ 2 is the Vandermonde determinant defined as: To evaluate the Riemannian unitary volume C R = U † dU , we use our previous technique [18] of evaluating a Gaussian integral in both rectangular and polar coordinates over an infinite domain. The integral is: From Riemannian volume invariance, this Gaussian integral can be evaluated by unitary transformation to an equivalent real antisymmetric matrix, which gives a onedimensional real Gaussian integral in each coordinate: Evaluating this using matrix polar coordinates: This has the form of a Mehta integral used in random matrix theory [46], and one obtains: Hence, on comparing with Eq (4.9) we find that the Riemannian unitary volume is given by: Since the matrices H used in Ref. [18] also have the class-D symmetry, we can relate this volume to the Euclidean volume C U for the class-D matrices derived previously [18]. As expected from Eq (4.4), the relation between these two factors is: We can now evaluate the invariant volume V = D dµ (ζ) of the classical domain using matrix polar coordinates. In this case we have to evaluate an integral of the form: (4.14) In order to perform this integral we consider the following change of variables: x j = ζ 2 j , hence dζ j = 1 2 x −1/2 j dx j . This allows the integral to have the form of a Selberg integral [46] 1 −1 (4.15) Using the result of C R given in Eq. (4.12) we obtain that the invariant volume V of the classical domain is:

C. Normalized basis
We now wish to prove that a resolution of unity for a suitably normalized Gaussian basisΛ N ζ is given by: (4.17) The Q-function basisΛ N ζ we shall use is a function of the stretched phase-space coordinates ζ on the classical domain of volume V, which vanishes at the domain boundaries. It is defined as: Here S ζ 2 is an even, positive scaling function, and N is a positive constant introduced for normalization purposes. This basis is positive definite, since the Gaussian operatorsΛ σ are positive definite [10,17], inside the classical domain. The function S must be invariant under unitary transformations, with a typical general form being: Using Grassmann variables, it is possible to prove that the normally ordered Gaussian operators can be brought into diagonal form by the unitary transformation of Eq (4.5) that diagonalizes the class-D hermitian matrix, and remain normally ordered in the new operator basis. This is shown in Appendix (A).
The trace normalization term for a Gaussian operator has an equal number of positive and negative eigenvalues. As a result, in terms of the eigenvalues of ζ and the transformed operatorsb = Uâ, the Gaussian operator is: (4.20) In order to obtain the resolution of unity we therefore must prove the following result, Here we have used the matrix polar coordinate Jacobian, Eq (4.6), the unitary volume of Eq. (4.12), and we introduce a radial integral:

D. Radial integral
To complete the proof of the resolution of unity, and obtain the normalizing factor N , we now focus on the radial part. This corresponds to the integral over the eigenvalues ζ. For simplicity in evaluating integrals we will take s = 0 in the normalizing function S ζ 2 . The normally ordered Gaussian operators can be expressed as: Using the result of Eq (4.23) and Eq (4.22), we can write the radial integral as: The value of the integralÎ ζ will depend on the normalization function. Expressed in terms of the eigenvalues ζ, this is: which is an even function of the eigenvalues.
The second integral of Eq. (4.24) has terms proportional to ζ j , which is an odd function of ζ j , while all the other terms are an even function of ζ j . Hence from the parity of these functions, the odd integrals vanish, soÎ ζ is proportional to a unit operator.
The radial integral which must be evaluated is: Next we perform the following change of variables: This allows the integral of Eq. (4.26) to be written in the form: which is another modified Selberg integral [46]. Therefore, the radial integral of Eq. (4.26) is: (4.27) Using the value of I ζ given in Eq. (4.27), together with the unitary integral, we obtain a general resolution of unity with a normalization constant N of the form: (4. 28) In the limit of k = 0, this normalization constant is just the phase-space volume V of the real classical domain, divided by the number of many-body quantum states, 2 M . This follows, since the requirement for the resolution of unity, in the limit of k, s → 0, is: The integral on the right hand side is the volume V of the real classical domain defined in Eq. (4.16). Hence in the limit of a uniform normalization, we simply have lim k,s→0 A Monte-Carlo integration was carried out to verify this result, by generating 10 6 random antisymmetric matrices X with −1 < X ij < 1, and testing positivity from the eigenvalues. This was in good agreement with Eqs. These results imply a simple physical interpretation of the resolution of unity. We note that for an arbitary operatorÔ, one immediately obtains from the k → 0 limit of the resolution, Eq (4.17), that: Hence, the average overlap of any Fermi operator with an orthonormal basis element equals its average overlap with a unit trace Gaussian operator.

V. FERMIONIC Q-FUNCTION AND OBSERVABLES
Following the method of Section II, the fermionic Qfunction is now defined inside the volume V in terms of the normalized Gaussian basis, as: This is the Hilbert-Schmidt inner product of the density operatorρ with the normalized Gaussian basis,Λ N ζ , defined in Eq. (4.18), and hence satisfies our first requirement.
The second requirement was that the Q-function should be a probability distribution which is normalized to unity. The density matrixρ, is hermitian and positive definite. The general normally ordered Gaussian operators defined in Eq. (4.18), are hermitian since the matrix ζ is hermitian, and they are also positivedefinite [10,17] within the classical domain. Thus, the fermionic Q-function is a non-negative probability distribution, and it is normalized from the resolution of unity.
As an example of a Q-function, the operatorÎ/2 M is the infinite temperature density matrix with unit trace. This has a constant overlap with any unit trace Gaussian of 2 −M , so that in the limit of k → 0 one obtains Q = 1/V, as expected for a uniform, normalized probability.

A. Differential identities of the fermionic basis
The third required property for a Q-function is that observables are moments of the distribution. In order to obtain this property, we will use the known differential identities given in Eq (3.10), which correspond to the action of the extended creation and annihilation operators on the Gaussian basis [10,17]. In particular, we will make use of the identity, â :â †Λ : = −σΛ −σ ∂Λ ∂σ σ.
Here the nested ordering {: . . . :} is defined as: The derivative of the fermionic Gaussian basisΛ N with respect to σ is obtained by using the product rule as well as this identity. Hence, the action of the ladder operators on the Gaussian basisΛ N for this particular nested ordering is given by: The extended creation and annihilation operators can be expressed in normal or antinormal form. We will consider the antinormal form of the observables, which are expressed as: Using the resolution of unity given in Eq. (4.17) we can write the observables as: Using the expression for the observables given in Eq. (5.6) as well as the differential identity of Eq. (5.4) and the definition of the Q-function given in Eq. (5.1), we obtain the following expression for the observables: Next, substituting the definition of the Q-function in this expression gives: On integrating the third term by parts, the boundary term vanishes due to the scaling factor S, for the corresponding bounded classical domain. We now consider the explicit form of the normalization function S and take the limit of k, s → 0 to simplify the result. The matrix derivatives are carried out on making use of Eq (3.11), which includes the class-D symmetry property. Hence, after changing to the stretched variables, the normally ordered Greens functions are given in a simple form by: where C M = 2M − 1/2 is a constant that depends on the dimensionality. This process can be iterated to obtain higher-order correlations.
In the infinite temperature case of Q = 1/V, the odd integral on the right-hand side vanishes, which means that â † iâ i = â iâ † i = 1/2. This is an expected result due to particle-hole symmetry.

VI. GAUSSIAN DENSITY OPERATORS
In this section we give explicit expressions of fermionic Q-functions in the case of Gaussian density operators. This includes the common cases of a non-interacting thermal gas and a BCS state. More generally, since the Gaussian operators are a complete basis for any density matrix, this allows any Q-function to be calculated if the corresponding fermionic P-function [10] is known. For simplicity, we will consider the limiting normalization with k, s → 0, unless stated otherwise.

A. Gaussian operator inner product
Consider two normalized Gaussian operators,Λ(σ) and Λ(σ ′ ). If we considerΛ(σ) as the Q-function basis, and Λ(σ ′ ) as a physical density matrix, then evaluating the Q-function reduces to evaluating their standard Hilbert-Schmidt inner product: In Appendix B, this is carried out following methods previously used to calculate linear entropy [34], together with Grassmann coherent states. This gives the result that: which takes a simpler form in terms of the stretched variables: If the anomalous terms vanish, so that m = m ′ = 0, then one obtains: 4) and this reduces to a previously obtained result [34] of: In the case of a general coordinate ζ, it is instructive to consider the overlap of a density matrix with itself. If we are on the limits of the classical domain, then I − ζ 2 = 0.
This implies that that all the eigenvalues have their extremal values of ζ j = ±1, which can only occur if there is a unitary transformation that sends every mode into either a ground or excited state, which is a pure state. As a consequence, we expect that F = 1, since F is just proportional to the linear entropy. This is exactly the result given by the inner-product expression, Eq (6.3).
Thus we see that the classical domain has a clear physical interpretation. It is a domain whose boundary is the Gaussian pure states, with the infinite temperature thermal state of highest entropy at the center.

B. Gaussian density operator Q-function
For any Gaussian density operator, like a thermal or BCS state, such that: the fermionic Q-function in the limit of k → 0 has the simple form: We have already given the simplest example of an infinite temperature state with ζ ′ = 0. The corresponding Qfunction is simply the uniform Q-function, Q = 1/V. More generally, since any density matrix can be represented [10] using a positive fermionic P-function P ζ ′ in the stretched coordinates, then any density matrix has a Q-function given by: A special behavior occurs in the case of a pure state, for which P ζ = δ ζ − ζ p , such that ζ 2 p = I. Let us call the corresponding pure state Q-function Q p . Given that it is a pure state ζ p , we may wish to evaluate the Qfunction at ζ = ζ p , which physically should be large from overlap arguments. The result can be found using the Q-function definition of Eq (5.1), together with the fact that a pure state density matrix is a projection operator. Alternatively, one can use the Gaussian state result of Eq (6.7). Either equation gives the same result: the Qfunction probability is increased by a factor of 2 M above its value in the highest entropy case: However, the Q-function is normalized to unity within the domain V. Therefore, the volume occupied by this high probability peak for a pure state must be relatively small at large M , of order 2 −M V.
In the opposite extreme, we can always find an antipodean state. This is exactly on the opposite side of the classical domain boundary, at ζ a = −ζ p . We note that ζ a ζ p = −ζ 2 p = −I. Hence, Q p ζ a = 0 for this pure state. That is, the Q-function of a pure state always vanishes at the point antipodean to the original pure state. Physically, this is caused by the fact that in the basis for which the pure state is a number state, the antipodean state has every eigenvalue reflected, and therefore it is orthogonal.
There are many other orthogonal states to a pure state. In fact, there are 2 M −1 orthogonal, pure boundary states ζ o of this type, in which one or more eigenvalues changes sign, leading to a vanishing Q-function with Q p ζ o = 0.
We note that by changing the values of k, s, the detailed shape of the Q-function can be modified. For example, if s ≫ 0 then the distribution is always concentrated close to the origin, and the high temperature state is a classical Gaussian similar to Eq (4.8). If s ≪ 0, the distribution is concentrated at the boundaries, and gives a distribution over pure states.
= ñ +â †â (2n − 1) . (6.10) Here we consider n real and n ∈ (0, 1); n gives the number of particles, whileñ = 1 − n gives the number of holes in the thermal state that corresponds toΛ 1 . As previously, it is useful to define a more symmetric variable, ζ = 1 − 2n ∈ [−1, 1]. Following the definition of the fermionic basis given in Eq. (4.18): where S 1 (n) is the corresponding normalization function of Eq. (4.19) for M = 1. If s = 0, the normalization function has a particle-hole symmetry, since it is an even function ofn, given by: For the single mode case the resolution of unity for the phase-space variable n is given by: Here we give an explicit proof for resolution of unity for the special case of one mode. Expanding the Gaussian, the integral over phase-space is: 14) The odd-parity integral terms overn all vanish, including the operator valued part. Hence, to prove a resolution of the identity, it is only necessary to show that: A resolution of unity therefore requires a normalization of: 16) in agreement with Eq (4.28). In the limit of k → 0, one has N = V/2 = 1.

D. Example of thermal states
We will now give an explicit expression for the singlemode Q-function described with this formalism. Since every single mode fermion state is a thermal state, we need only consider the thermal states, with the finitetemperature Fermi-Dirac mean occupation number of: The density matrix is itself one of the Gaussian operators: The explicit form of the fermionic Q-function for this case is given by taking the inner product of Eq. (6.18) and the normalized Gaussian stateΛ N (ζ). The symmetrized thermal occupation is ζ th = 1−2n th = 2ñ th −1. Therefore, in the k, s → 0 limit, the Q-function is given by: Here we have used the result for the trace of two normally ordered normalized Gaussian operators [34]. Since −1 < ζ, ζ th < 1, the Q function is clearly positive, and the result agrees with the general expression for a multi-mode Gaussian Q-function, Eq (6.7).

E. Q-function and observables
Using the expression for the observables given in Eq. (5.2) , we can write the corresponding differential identities as:ââ †Λ Hence, since the normalization factor is N = 1, the observables for the single mode case are: On performing the derivatives, taking the limit of k → 0 and simplifying terms, we obtain: Using the expression of the Q-function given in Eq. (6.19) and the expression for the observables given in Eq. (6.22), we evaluate the observable, which is This result corresponds to the expected thermal average value.

VII. SUMMARY
In summary, we have introduced a probabilistic fermionic Q-function for an arbitrary many body fermionic system. The fermionic Q-function is physically interpreted as a suitably normalized overlap of the density operator and the normally ordered Gaussian operators.
We have also used three important properties of these Gaussian operators, which are their positivity inside a bounded domain, their resolution of unity and their differential identities. As a result, we have shown that the Q-function is a continuous probability distribution, exists for any quantum density matrix and has observables which are moments of the distribution.
The fermionic Q-function derived here is a general probabilistic approach to fermion physics. It uses as a complete basis the Gaussian operators, which have been applied in the study of many-body fermionic systems. This is in contrast to previous fermionic Q-functions. For example, a Q-function defined in terms of anticommuting Grassmann variables is neither a probability nor even a real number.
There is an elegant physical interpretation of the resulting phase-space domain. These Gaussian states are physical states which have the highest possible entropy at the domain centre, and the lowest possible entropy at the boundaries. This is because the classical domain boundary is the set of Gaussian pure states, which includes the BCS states. On the other hand, its centre at ζ = 0 is the infinite temperature thermal state of maximal entropy.
Thus, as a hot fermionic system is cooled, we expect its Q-function to be initially uniform on the classical domain, while gravitating towards the boundary as it cools towards states of greater purity.

Gaussian series expansion
We wish to prove thatÎ Λ , and henceΛ u , is equal to a normally-ordered diagonal Gaussian form. To proceed, Λ u inside the integral is expanded as a series of even order normally-ordered polynomials in the Fermi operators, so that:Λ Using an explicit form of the matrix µ and recalling thatâ = â,â † , we can normally order the lowest order quadratic form P (2) ofΛ (u) as: In the last step there is a sign change as well as a transposition, owing to the definition of fermionic normalordering, so that :bâ † := −â †b . From the Grassmann coherent state eigenvalue properties and anticommutation relations, we obtain: and applying this to all terms: where Σ is the transposition matrix defined in the main text. For higher order polynomials, an analogous procedure occurs. There will be an even number of sign changes for both Fermi operator re-orderings and for Grassmann variable re-orderings. Hence:

Unitary transformations
So far we have expressed the un-normalized Gaussian operator term in terms of Grassmann variables. Next, we are going to consider the class-D unitary transformations that diagonalize the matrix µ under consideration. This is a 2M × 2M unitary transformation of the Fermi operators: Note that if we diagonalize the matrix µ we also diagonalize the matrix σ. As pointed out by Altland and Zirnbauer [19], U must satisfy U = ΣU * Σ with: This property implies that: which guarantees that the transformation, when applied to the Fermi operators, preserves operator commutators [38]. We can apply the same transformation to the Grassmann variables, by choosing that From the first expression, one obtains: while on conjugating the second expression: α ′ * = u 12 * β + u 11 * α = u 21 β + u 22 α = α ′ .
This means that the transformation is a consistent transformation on Grassmann variables that preserves conjugation, so that if |α ′ ′ is the new coherent state in the new basis, then: Therefore, the Grassmann quadratic in the exponential is also diagonalized by the unitary transformation, so that: It follows that one can extend this argument to each even order polynomial form in the exponential, so that: exp −γ T Σµγ/2 = exp −γ ′T Σλγ ′ /2 . (A20)

Gaussian in operator form
We can now evaluate the Grassmann integral, noting that the integration measure is invariant under a unitary transformation, since it has unit Jacobian: We can now use the eigenvalue properties of the Grassmann coherent states in the transformed basis to return back to the operator form in the diagonal basis, so that, on making use of the Grassmann coherent state expansion of the identity, one obtains: This means that we have the required result. That is, we can diagonalize the class-D covariance matrix for the normally-ordered Gaussian operators using the same transformations used in [18], while leaving the underlying normally-ordered form invariant, just as one can in a non-normally-ordered case.
Appendix B: Inner product of two Gaussian operators Consider two normalized Gaussian operators,Λ(σ) andΛ(σ ′ ). We wish to evaluate their standard Hilbert-Schmidt inner product: F σ, σ ′ = Tr Λ (σ)Λ(σ ′ ) . (B1) We will start be evaluating the inner product of the two corresponding un-normalized operators,Λ u µ and Λ u (µ ′ ), where: We now use the known result for a trace, and the expansion of the identity operator in Grassmann variables, so that:

Sign reversal variable change
It is convenient at this stage to make a variable change ofᾱ = −ᾱ ′ , noting that in Grassmann calculus, the conjugate variable can be regarded as an independent complex variable. This causes a sign change in the integration of (−1) M . As a result, we can introduce a new variable, and its conjugate, which is: We note that d 2M αd 2M β = (−1) M d 2M αd 2Mβ = d 4M γ, since the sign change and the Grassmann reordering in the integration variables will cancel each other. As a result, we can write that where Γ is a 4M × 4M matrix, such that:

Grassmann integration result
This is a standard Grassmann gaussian integral, except in a space of twice the previous dimension. Hence, the inner product is a Pfaffian of the antisymmetrized form of Γ, given as determinant by: Now, changing variables to the covariance matrix, one obtains:Ī so that in terms of these variables, Hence, multiplying by the normalizing factors of det iσ , and noting that (−1) 2M = 1, one obtains: This can also be written in terms of the stretched variables, ζ =Ī − 2σ. We note that: