On the partial transpose of fermionic Gaussian states

We consider Gaussian states of fermionic systems and study the action of the partial transposition on the density matrix. It is shown that, with a suitable choice of basis, these states are transformed into a linear combination of two Gaussian operators that are uniquely defined in terms of the covariance matrix of the original state. In case of a reflection symmetric geometry, this result can be used to efficiently calculate a lower bound for a well-known entanglement measure, the logarithmic negativity. Furthermore, exact expressions can be derived for traces involving integer powers of the partial transpose. The method can also be applied to the quantum Ising chain and the results show perfect agreement with the predictions of conformal field theory.


Introduction
Entanglement plays a key role in the study of quantum many-body systems [1,2]. Considering a pure state of a composite system, a simple measure of the entanglement between two complementary parts is given by the von Neumann (or entanglement) entropy. Particularly interesting is the case of pure ground states where, in great generality, an area law for the entanglement emerges [3]. The most well established exceptions are one-dimensional quantum chains at criticality, where the entanglement entropy shows a universal logarithmic scaling [4] which can be fully understood with the help of conformal field theory (CFT) [5]. The predictions of CFT have been confirmed on a variety of lattice models, among which a distinguished role is played by free-particle Hamiltonians. The ground states of these systems are given by bosonic/fermionic Gaussian states where the full entanglement spectrum is easily accessible [6].
The characterization of entanglement for mixed states is, however, far less obvious since, in contrast to the pure-state scenario, there is no unique way of defining a wellbehaved measure. Among the numerous proposals for entanglement measures [7], a large family is based on a convex-roof extension of the von Neumann entropy. The drawback of these constructions is that they are essentially uncomputable already for systems of relatively small size. A viable alternative is based on an entirely different approach, making use of a special property of the partial transposition. Namely, the spectrum of the partial transpose of a density matrix may contain negative eigenvalues, only if the state is entangled [8,9]. In turn, a measure called logarithmic negativity [10] can be introduced, which quantifies how much the partial transpose of a state fails to be positive, and can be shown to fulfill all the requirements of an entanglement measure [11].
Although being a computable measure, the evaluation of logarithmic negativity might still pose a significant challenge in extended quantum systems. A notable exception is the case of bosonic systems, where the effect of partial transposition is equivalent to a partial time-reversal of the momenta in the corresponding subsystem [12]. Furthermore, the partial transpose of bosonic Gaussian states remains to be Gaussian and, in turn, one has a simple formula to compute the logarithmic negativity via the covariance matrix [13]. Remarkably, the analogue statement does not hold for fermionic Gaussian states.
The early studies of logarithmic negativity in lattice systems were conducted for the harmonic oscillator chain [13,14,15,16,17,18] using the covariance matrix technique, and for spin chain models [19,20] via density matrix renormalization group calculations. In contrast, exact analytical results were found only for a few simple spin models [21,22]. A renewed interest in the problem was triggered recently, after a systematic approach within CFT was introduced [23]. This method could be applied to calculate the entanglement negativity for various geometries in ground [24,25] or thermal states [26,27] of one-dimensional systems, as well as in out-of-equilibrium situations [27,28,29,30].
Even though the predictions of CFT can be routinely tested on harmonic chains, calculating the logarithmic negativity in fermionic or spin systems remains to be more difficult. Recent studies employed a tensor-network representation of the partial transposition to calculate entanglement negativity for the transverse Ising chain [31]. Alternatively, Monte Carlo techniques were applied to calculate higher moments of the partial transpose [32,33]. However, even for the simplest case of fermionic Gaussian states, a method which could compete with the computational simplicity of the bosonic case has so far been unknown.
Here we show that, with a suitable choice of basis, the partial transpose of fermionic Gaussian states can be cast in a particularly simple form. Namely, it can be written as the linear combination of only two Gaussian operators, uniquely defined by corresponding covariance matrices which can be found explicitly. Under further assumption of a reflection symmetric geometry, this construction can be used to calculate a lower bound for the logarithmic negativity via the covariance matrix spectrum. For critical systems, the scaling behaviour of this bound shows remarkable similarities to that of the entanglement negativity. Furthermore, the higher moments of the partial transpose can be exactly evaluated through simple trace formulas, providing a way to test the universal CFT predictions on fermionic Gaussian states with minimal computational costs.
In Section 2 we define fermionic Gaussian states and introduce the specific models in consideration. The partial transposition transformation for fermions is discussed in Sec. 3, focusing on a particular choice of basis which leads directly to our main result. Sec. 4 is devoted to the construction of a lower bound for the logarithmic negativity, and its numerical investigation for the quantum Ising chain. Trace formulas for integer powers of the partial transpose of the reduced density matrix are presented in Sec. 5. The paper concludes in Sec. 6 with a short discussion of the results and their possible extensions. Various details of analytical calculations are included in three Appendices.

Model and definitions
We consider quantum systems associated to free-fermion Hamiltonians where the matrices A and B are Hermitian and antisymmetric, respectively. The fermionic creation/annihilation operators, c † m and c m , satisfy the canonical anticommutation relations c † m , c n = δ mn . For our purposes, it will sometimes be more convenient to work with Majorana fermions defined as satisfying the relations {a k , a l } = 2δ kl . In terms of Majorana operators, the Hamiltonian of Eq. (1) with real A and B can be rewritten as H = i 2N m,n=1 T m,n a m a n , where T 2m,2n−1 =−T 2n−1,2m = 1 4 (A m,n +B m,n ) and T 2m,2n =T 2m−1,2n−1 =0. The product of all Majorana operators define the parity operator, P = i N 2N n=1 a n , which plays an important role in fermionic systems. According to the parity superselection rule, only density matrices that commute with P correspond to physical states [34,35,36].
The states we are going to study in this paper are the so-called Gaussian states. These describe the ground and Gibbs states of quadratic Hamiltonians, and play a prominent role in quantum information theory [37,38,39]. A state ρ is Gaussian if it can be written as where W is a purely imaginary antisymmetric matrix (with the possibility of |W kl | → ∞ allowed) and Z is the normalization factor. A Gaussian state can also be characterized uniquely by its covariance matrix Γ kl = [a k , a l ] /2 via Using the covariance matrix, one can express the expectation value of any Majorana monomial through the Wick expansion Tr(ρ a n 1 a n 2 . . . a n 2ℓ ) = π sgn (π) 2ℓ k=1 Γ n π(2k−1) ,n π(2k) , where all the indices are different, the sum runs over all pairings π, and sgn (π) denotes the sign of π ‡. Let us note, that similarly to Gaussian states, one can introduce general Gaussian operators which are also defined through Eq. (3), however, without requiring that the spectrum of W is real. The Wick expansion, i.e. Eq. (5), holds for these operators, as well. We will also study spin chain models that are related to free-fermion Hamiltonians through the Jordan-Wigner transformation [40] where σ α m (with α = x, y, z) denote the Pauli matrices acting on site m. The prototypical example is the transverse field Ising (TI) chain, where a chain of length N with open boundary conditions is considered. Applying the Jordan-Wigner transformation, the TI Hamiltonian (7) takes the form of Eq. (1) with ). (8) ‡ A pairing π over a set {1, 2, . . . , 2ℓ} is a permutation of the 2ℓ elements which satisfies π(2k − 1) < π(2k) and π(2k − 1) < π(2k + 1) for any 1 ≤ k ≤ ℓ. Any pairing can be decomposed into an M number of transpositions (simple exchange of only two indices), and the sign of the pairing is defined as sgn (π) = (−1) M .
Such a quadratic Hamiltonian can be diagonalized by a canonical transformation, and brought into the standard form The spectrum Λ k and the vectors φ k , ψ k in Eq. (9) follow from the eigenvalue equations With the full solution of the problem at hand, one can directly write down the covariance matrix for a Gibbs state, with inverse temperature β, of the open TI chain as with matrix elements given by In our studies we will also be interested in the periodic TI chain, given by a Hamiltonian as in Eq. (7) with the sum running up to L and boundary condition σ x L+1 = σ x 1 . Note that, for clear distinction, we will use L instead of N for the length of the periodic chain. It is well known, that its ground state is the same as for the fermionic model with matrices as in (8) but with antiperiodic boundary conditions. The solutions of the system (11) and (12) are plane waves, φ k (m) ∼ e ip k m and ψ k (m) ∼ e iθ k e ip k m , with the Bogoliubov angles and eigenvalues given by and the allowed values of the momenta are p k = (2k−1)π/L with k = −L/2+1, . . . , L/2. One can also work directly in the thermodynamic limit, L → ∞, where the momenta become continuous and the sum in Eq. (14) for the matrix elements g mn is replaced with an integral.

Partial transpose for free fermions
As discussed in the Introduction, the partial transposition plays an important role in quantum information theory. In the context of entanglement theory, it was first studied for qubit and qudit systems [8,9], but later also bosonic [12,41,42,43] and fermionic models [34,35,44] were investigated. An important result coming from these studies was that the partial transpose of a bosonic Gaussian state is again a Gaussian operator; this simplifies the calculation of the negativity [13,45]. The analogous result for fermionic Gaussian states does not hold, which can already be demonstrated by 2-site systems, see Appendix A. This section is devoted to the derivation of a weaker, but still useful, result for the fermionic case. After recalling the notion of the partial transpose for spin systems and the corresponding definition for fermions in Section 3.1, we show in Section 3.2 that the partial transpose of a Gaussian state (in a particular basis) can always be decomposed as a sum of two Gaussian operators. This decomposition lies at the heart of all the results in the further sections. Consider a general tripartition of a chain of qubits into disjoint sets A 1 , A 2 and B, e.g. as in Fig. 1. Let ρ denote the density matrix of the whole composite system. Defining A = A 1 ∪ A 2 , the reduced density matrix (RDM) of subsystem A is given by ρ A = Tr B ρ. The partial transpose of the RDM, ρ T 2 A , with respect to the subsystem A 2 is defined by its matrix elements as where {|e } denote complete bases on the Hilbert spaces H 1 and H 2 pertaining to the subsets A 1 and A 2 . The definition of ρ T 2 A is basis dependent. However, one can easily characterize the set of transpositions on the operators acting on H 2 as those non-degenerate linear transformations that satisfy for any two operators M 1 and M 2 acting on H 2 . Since any two partial transpositions can be connected by a unitary conjugation, the eigenvalues of ρ T 2 A are independent of the choice of basis. Moreover, it was shown that the partial transpose of the density matrix can only have negative eigenvalues if the corresponding state is entangled [8,9].
In a similar way, one can define the partial transpose for fermionic states. Consider a tripartition of a system with N fermionic modes, e.g. as in Fig. 1 . Let {m 1 , m 2 , . . . , m 2k } and {n 1 , n 2 , . . . , n 2ℓ } denote the indices of the Majorana operators belonging to the subsystems A 1 and A 2 , respectively. Let us introduce the notation a 0 A general fermionic state on A = A 1 ∪ A 2 can be written as where the variables κ = (κ 1 , . . . κ 2k ) and τ = (τ 1 , . . . , τ 2ℓ ) in the summation run over all bit-strings of length 2k and 2ℓ, respectively. Note that since physical fermionic states must commute with the parity operator, as discussed in Section 2, one has w κ,τ =0 when 2k i=1 κ i + 2ℓ j=1 τ j is odd. The partial transpose of ρ A is simply a transformation that leaves the A 1 Majorana modes invariant and acts as a transposition on the operators built up from modes of where R satisfies Eq. (17). Since also in the fermionic case all the transpositions are connected by a unitary conjugation, the eigenvalues of ρ T 2 A will be independent of which R we choose. It will be useful to consider the following particular transposition which is defined by where |τ | = 2ℓ i=1 τ i . In other words, a Majorana monomial is mapped by R to itself if it is of length 4n or 4n + 3, and otherwise it is multiplied by a −1 sign. Note that a very similar transposition in fermionic systems was already considered in Ref. [46]. Although the definition of entanglement in fermionic systems is somewhat different from the case of spin systems, it has been proven that ρ T 2 A can only have negative eigenvalues if ρ A is entangled also in the fermionic case [34,35].
Finally, let us shortly discuss the connection between reduced density matrices of fermionic and spin models that are connected by the Jordan-Wigner transformation. As this transformation is non-local, it has been shown that the reduced density matrices corresponding to a region A 1 ∪ A 2 in a spin chain model and its fermionic counterpart are usually not equivalent (not isospectral), unless A 1 and A 2 are adjacent intervals, as depicted in Fig. 1(b) [47,48]. Since the same holds also for the transposed density matrices, when treating spin models we will only consider the adjacent interval geometry. For the case of fermionic systems, our results are valid for arbitrary geometries.

The Gaussian case
Consider a Gaussian state ρ A on the system A = A 1 ∪ A 2 with a covariance matrix where Γ 11 and Γ 22 denote the reduced covariance matrices of subsystems A 1 and A 2 , respectively; while Γ 12 and Γ 21 contain the expectation values of mixed quadratic terms. Let P 2 be the parity operator on subsystem A 2 , and define the operators . By definition, we have that ρ A = ρ + +ρ − . Using the notation of Section 3.1, ρ + and ρ − can be expanded as where the coefficients w κ,τ can be obtained from Γ A using the Wick rule, Eq. (5). By linearity of the partial transpose, ρ T 2 A = ρ T 2 + + ρ T 2 − follows, and ρ T 2 ± can be obtained using Eq. (20): Let us introduce the generalized Gaussian operators O + and O − , with covariance matrices and consider the Majorana monomial expansion of these operators, Since O + and O − are Gaussian operators, one can again obtain o ± κ,τ from Γ ± using Eq. (5). Connecting the Wick-expansion form of w κ,τ with that of o ± κ,τ , using the relation between Γ A and Γ ± , one can deduce that Comparing this with Eq. (23) it immediately follows that ρ T 2 . Thus, we obtain the decomposition which is, from a conceptual point of view, the main result of the paper.

Partial transpose and logarithmic negativity
In the previous section, we have shown that the partial transpose of a Gaussian RDM can be written as a linear combination of only two Gaussian operators, which is the simplest possible form for a non-Gaussian operator. However, since O + and O − do not commute in general, Eq. (27) can not be rewritten for the eigenvalues, and thus one does not have a direct access to the full spectrum of ρ T 2 A . Nevertheless, there are a number of important properties which can be deduced by a direct investigation of the covariance matrices Γ ± . A particularly interesting quantity that we will study is the logarithmic negativity [10], which can be used as a measure of entanglement.
In the following we thoroughly investigate three special cases. First, we consider the partial transpose for bipartite pure states which, although the results being well-known, turns out to be very instructive in understanding the implications of the decomposition in Eq. (27). We proceed with the study of thermal mixed states in a reflection symmetric bipartite geometry, which allows us to define and calculate a lower bound for the logarithmic negativity. Finally, we report our findings for a genuine tripartite geometry. In each of the following subsections, the validity of formula (27) was checked against exact numerical calculations for the TI chain, Eq. (7), with a small number of spins.

Bipartite pure states
We first consider the simplest case with B = ∅ and a pure state on A = A 1 ∪ A 2 given by ρ = |φ φ|. Since for any pure Gaussian fermionic state Γ 2 = ½ is satisfied, it follows that [Γ + , Γ − ]=0, which implies that O + and O − commute as well. Furthermore, since their eigenvalues are connected by a complex conjugation, the spectrum of ρ T 2 is simply given as the sum of the real and imaginary parts of the O + eigenvalues. Now, since O + is also Gaussian, its spectrum is uniquely determined through the eigenvalues of Γ + . To obtain them, we first assume without loss of generality that |A 1 | ≤ |A 2 |. Furthermore, for notational simplicity we also assume that |A| is even, however, the results for the odd case follow trivially. Let ±µ k denote the eigenvalues of the reduced covariance matrix Γ 11 with k = 1, . . . , |A 1 | and µ k ≥ 0. As shown in Appendix B, the eigenvalues ±ν ± k of Γ + can then be given as The canonical diagonalised form of the Gaussian operator O + reads where b ± j are Majorana operators obtained from a j via the orthogonal transformation which diagonalizes Γ + . The eigenvalues of O + can be obtained according to the following rules. First, we shall consider the various combinations of the conjugate eigenvalue pairs for each k = 1, . . . , |A 1 | as with σ k = ± and σ ′ k = ±. Using Eq. (28) this yields where we used the property ν + k ν − k = 1. The nonzero eigenvalues of O + can then be written down as where σ and σ ′ are the signature arrays corresponding to the eigenvalue. Note that the additional ν ± k = 1 in Eq. (28) lead to further eigenvalues of O + that are all equal to zero.
The products in Eq. (32) are either real or purely imaginary and the eigenvalues of ρ T 2 thus follow as Re Ω σ σ ′ or Im Ω σ σ ′ , respectively. It is instructive, however, to derive the same spectrum using the Schmidt decomposition of |φ . Dividing the Hilbert space with the RDM eigenvalues λ i of subsystem A 1 . Clearly, the state is supported on a smaller Hilbert space H 1 ⊗ H 1 and is invariant under the action of a flip operation defined by |φ 1 i |φ 2 j → |φ 1 j |φ 2 i . The partial transpose of ρ, commutes also with the flip operator. Furthermore, it is easy to check that the eigenvalues and vectors are where we introduced the notation Note that all the positive (negative) eigenvalues correspond to even (odd) eigenvectors with respect to the flip operation. Moreover, since ρ is Gaussian, one can immediately write down the products of eigenvalues as where the signature σ has again components σ k = ±. Comparing Eq. (37) to Eq. (32), one indeed recognizes Ω σ σ ′ up to the factors of i.
Owing to the simple product structure of the eigenvalues in Eq. (37) and, in particular, to the fact that all the negative eigenvalues are located in the odd subspace, one has (38) § When considering a tensor product structure, we always refer to a partition of the spin chain system constructed through the Jordan-Wigner transformation.
Thus, the pure-state logarithmic negativity is given by the simple formula Considering also the expression of the Rényi entropy for fermionic Gaussian states, the well-known equality E = S 1/2 for pure states can be confirmed directly.

Thermal states in a bipartite geometry
The simplicity of the pure-state scenario relies essentially on the property of the Schmidt decomposition, which is automatically symmetric under the flip operation defined previously. Due to this, the partial transposed state is block diagonal wrt the splitting into even and odd subspaces. Such a structure is missing for general bipartite mixed states, unless the system has a flip-type symmetry a priori. In this respect, a natural scenario would be to consider intervals of equal length |A 1 | = |A 2 | = N/2 and states that are reflection symmetric. To analyse such a situation, we shall consider Gibbs states of the open TI chain with a symmetric bipartitioning, these being the simplest mixed Gaussian states where we hope to get further insight into the structure of the partial transpose.
As the spectrum of Γ + is invariant with respect to a sign change and complex conjugation, hence the eigenvalues can be collected into two families of quadruplets (41) where in family (I) we choose Re z k > 0 and Im z k > 0, whereas u k > v k in family (II). Note that the eigenvalues in the second family are purely imaginary and thus their complex conjugate are automatically contained in the spectrum of a skew-symmetric matrix, hence u k = v k . Although one could, in general, have an arbitrary number of type (II) quadruplets, from the numerics we observe that in the Ising case they are either absent or a single one appears. Moreover, this only happens in the symmetry-broken phase, i.e., when h < 1 in Eq. (7). Analogously to the pure case in Eq. (30), we first assign the factors within each quadruplet, and the eigenvalues Ω σ σ ′ of O + are again given in the factorized form of Eq. (32). Although the spectrum of the operator O − is identical to that of O + , they do not commute in general and thus one has no direct access to the eigenvalues of ρ T 2 . Nevertheless, the information about the even/odd parity of the eigenvectors is retained. The reflection operator R, which defines the even/odd subspaces in our case, acts on the spin operators as Rσ α n R † = σ α N −n (with α = x, y, z), implying the action Rc † n R † = P c † N −n on the creation operators, where P is the parity operator. Using this, it follows that the sign factor associated to an eigenvector reads which can also be verified by considering the pure state limit.
With the knowledge of S σ σ ′ , we are now able to carry out signed traces of the form where Tr e/o denotes the trace taken over the even/odd subspace. Note that the terms in the sum of Eq. (44) Finally, we define the quantity which clearly gives a lower bound for the logarithmic negativity, E ≥ E o , with strict equality if all the negative eigenvalues reside in the odd sector and their number is equal to the dimension of that subspace. This is true for pure states and one expects it to be valid for thermal states in a finite regime of temperatures above the ground state. For general bipartite states, however, the dimension of the negative subspace can be much larger [49,50].
To test the bound E o against the exact value of the logarithmic negativity, we considered small TI chains with N ≤ 10 and obtained E via exact diagonalisation of ρ T 2 . This is shown on Fig. 2, as a function of the temperature (left) as well as of the magnetic field (right). We find that, for low enough temperatures, E o indeed exactly coincides with E. For larger temperatures, however, some of the negative eigenvalues in the odd sector become positive or, vice versa, even eigenvectors could attain negative eigenvalues, with both of these processes increasing the difference E − E o . Nevertheless, for the small system sizes considered, it appears that E o gives a very good approximation up to temperatures T ≈ 1/N, above which it starts to deviate significantly.

Ground states in a tripartite geometry
So far we have only considered bipartite geometries. Another interesting setting we are able to deal with is a pure state with the symmetric tripartite geometry, depicted on Fig. 1(b). In this case, the reduced state ρ A after tracing out the sites of B is a mixed Gaussian state, associated to the reduced covariance matrix Γ A , with indices running over sites in A [51]. The logarithmic negativity for the tripartite case can be obtained with CFT methods [23]. For two intervals of the same size ℓ embedded in a system of length L with periodic boundary conditions, the calculation yields [24] E(ℓ, L) = c 4 ln L π with the central charge c and a non-universal constant. However, subtracting the value at ℓ = L/4 one obtains a universal scaling function where z = ℓ/L and we have set c = 1/2 corresponding to the TI chain. The formula was tested using tensor network methods for the calculation of the partial transpose, and a very good agreement was found [31]. It is interesting to check the behaviour of the lower bound, defined in Eq. (46), for the geometry at hand. In fact, since reflection symmetry is fulfilled, all the arguments of the previous section, leading to Eq. (45), apply and E o is given by the same formula in terms of the eigenvalues of Γ + . However, there is an important difference compared to the bipartite thermal case, which is apparent from the numerical investigation of small systems. Namely, the number of negative eigenvalues of ρ T 2 A is always less then the dimension of the odd subspace and, moreover, some of the corresponding eigenvectors are even. Thus, by tracing out the sites of B, the partial transpose ρ T 2 A cannot be smoothly deformed from the pure-state case and, consequently, one does not have a finite regime of parameters where the bound given by E o is tight.
In spite of the above findings, E o shows a very interesting behaviour, which is demonstrated on Fig. 3. First of all, it shows a clear signature of the phase transition at h = 1, which can be seen when plotting E o (L/4, L) against h on the left of Fig. 3. Furthermore, defining the quantity ǫ o analogously to Eq. (48), one finds an excellent data collapse when plotted against the variable z = ℓ/L, see right of Fig. 3. The scaling function is found to be given by Although the functional form of ǫ o (z) was found by trial, one has an excellent match with the data without any fitting parameters involved. Interestingly, the prefactor of the logarithm is exactly the half of ǫ(z) in Eq. (48), however, the argument is modified as well. We also performed a calculation directly in the thermodynamic limit L → ∞, and found E o (ℓ) = 1/16 ln ℓ + const., which is perfectly consistent with the above findings. Furthermore, one could also consider the simple fermionic hopping chain (or XX chain in spin language), defined by B mn = 0 and A mn the same as in Eq. (8). In complete analogy with the result for the bipartite entanglement [52], one finds E XX o (2ℓ) = 2E T I o (ℓ) for h = 0, and thus a doubled prefactor 1/8 with respect to the TI chain. Therefore, even though E o does not approximate E well, it shows exactly the same universal behaviour, suggesting it as an entanglement indicator which is extremely simple to calculate.

Trace formulas
Our result for the partial transpose, Eq. (27), can be further tested by looking at traces of integer powers of the partial transpose ρ T 2 A which can also be carried out within CFT. Since the identity Tr(ρ T 2 A ) 2 = Tr ρ 2 A holds for any density matrix, the simplest nontrivial quantity to check is the trace of the third power. For the geometry of the previous section, one finds the CFT result [24] R 3 (ℓ, L) = ln Tr(ρ T 2 A ) 3 = − 1 9 ln L 3 π 3 sin 2 πℓ L sin 2πℓ L + const. (50) Similarly to Eq. (48), a universal scaling function can be defined as [24] which was already tested numerically for the critical TI chain [31,32]. On the other hand, one could also consider two adjacent intervals of equal length ℓ, embedded in an infinite chain which is thermalized at inverse temperature β. Applying a simple conformal transformation, the corresponding CFT formula follows as We now show how the above traces can be calculated with the covariance matrix formalism. Expanding the third power of ρ T 2 A in Eq. (27) and taking the trace one arrives at where we have used that both of the traces on the right hand side are real. In order to evaluate them, one has to invoke the determinant formulas for the trace of products of Gaussian operators, which have already been considered in different contexts [47,53].
The main steps of this calculation are summarized in Appendix C. In turn, one finds where the sign of the first term depends on the spectrum of Γ + , see Eq. (41). Namely, the + sign applies only if the quadruplet with purely imaginary eigenvalues appears (see Appendix C for a more detailed discussion). Note that similar sign ambiguities also appeared in [47]. One should also remark, that traces of higher powers can be handled in a very similar way, however, the formulas become rather lengthy. The trace formula (54) can now be compared to the CFT predictions in (51) and (52) by inserting the corresponding covariance matrices Γ ± and evaluating the determinants. This is shown in Fig. 4 for the ground (left) and thermal states (right), respectively. The perfect agreement of the curves provides a highly nontrivial check of the CFT results.

Discussion
In conclusion, we have shown that the partial transpose of fermionic Gaussian states can be written as a linear combination of only two Gaussian operators, uniquely determined by associated covariance matrices Γ ± . In the presence of reflection symmetry, this particular form of the partial transpose allows us to carry out traces over the even/odd subspaces which, in turn, can be used to construct a lower bound to the logarithmic negativity. Furthermore, the trace of any integer power of ρ T 2 A can, in principle, be calculated as a sum of determinants, each of linear size 2|A|.
There are several open questions left for future research. We did not consider in detail entanglement detection questions, e.g., providing temperature bounds for separability of fermion or spin systems in thermal equilibrium. It would be instructive to compare such results obtained from the negativity lower bound E o with earlier studies [54,55,56].
Another natural extension of our work would be to consider non-adjacent intervals. For spin chains, however, it was shown that the spin RDM itself is already a linear combination of four fermionic RDMs [47]. Although our construction for the partial transpose could be carried over, it would further double the number of terms in the linear combination. Thus, the calculation of the traces for such a case is still realisable, but presumably more tedious.
It would also be interesting to see whether the lower bound E o could be attainable within the framework of CFT. This could lead to an analytical understanding of the scaling function for the critical tripartite case in Eq. (49) and could shed light to the origin of the prefactor. In fact, one is tempted to guess that this is equal to one-half of the corresponding prefactor of the logarithmic negativity in a general CFT. From the free-fermion point of view, our analysis clearly suggests that certain asymptotic relations between E o and E could hold in general. Finding a rigorous form of this relation would allow for a numerically feasible estimation of the entanglement negativity for fermionic systems.
The work of V.E. was supported by OTKA Grant No. NK100296; and Z.Z. acknowledges funding by the British Engineering and Physical Sciences Research Council (EPSRC).
Appendix A. The partial transpose of a 2-site RDM It is instructive to check how the method introduced in Section 3 works for the simplest case of two consecutive sites. The canonical form of the Gaussian RDM is given by where ν k ∈ [0, 1] and the Majorana operators b j are related to a j through an orthogonal transformation. For simplicity, we will consider only covariance matrices of the form of Eq. (13), and assume also reflection symmetry. These states are parametrized by a single angle θ beside the covariance matrix eigenvalues ν k . Using the Jordan-Wigner representation of a j and working in the usual spin basis, the most general form of the partial transpose is Note that, besides the diagonals, all matrix elements vanish and are thus not shown. It is straightforward to obtain the four eigenvalues λ 1,2 = 1 + ν 1 ν 2 ± (ν 1 + ν 2 ) 2 cos 2 2θ + (ν 1 − ν 2 ) 2 /4, Using the parity operator P 2 = σ z 2 , one can also construct ρ T 2 ± = (ρ T 2 ± P 2 ρ T 2 P 2 )/2 with matrix elements The eigenvalues of the operator ρ T 2 + + iρ T 2 − then read , Note that we have λ 3,4 = Re Ω 3,4 + Im Ω 3,4 which simply follows from the fact that ρ T 2 + and ρ T 2 − commute on the corresponding subspace, including the odd eigenvector. Unfortunately, this property does not generalize to symmetrically bipartitioned intervals with more than two spins.
We will now show that the Gaussian operator O + with covariance matrix Γ + has indeed eigenvalues given by Eq. (A.6). The covariance matrix Γ for the Gaussian state (A.1) and the associated Γ + have the form with the shorthand notation The four eigenvalues ±ν ± of Γ + can be computed with If the operator O + is Gaussian, its eigenvalues must have the form Finally, let us shortly discuss the non-Gaussian character of ρ T 2 . comparing the expectation values Tr(ρ T 2 a m a n ), where m, n=1, . . ., 4, with Tr(ρ T 2 a 1 a 2 a 3 a 4 ), one observes that the Wick expansion, Eq. (5), does not hold, unless c 2 = ν 1 ν 2 . This remains true whatever basis we choose for the partial transpose. Thus, the partial transpose of a fermionic Gaussian state is usually not a Gaussian operator.

(B.5)
Thus the spectrum of Γ ± follows from the eigenvalues of the squared matrix Rewriting in terms of the eigenvalues (ν ± k ) 2 and µ 2 k of G ± and G, respectively, one arrives at (B.11)

Appendix C. Determinant formulas
Let us consider the Gaussian operators O ± corresponding to the generalized covariance matrices Γ ± = tanh(W ± /2). With a denoting the vector of Majorana operators, we introduce O (W ± ) = exp aW ± a 4 , Z(W ± ) = Tr exp aW ± a 4 , (C.1) and thus O ± = O (W ± ) /Z(W ± ) Our aim is to calculate various traces of the form Tr(O m + O n − ) with some integers m and n. Following the lines of Ref. [47], we first introduce the notation We also note the simple fact that O m (W ± ) = O (mW ± ). Hence, the traces we consider can be written in the form They can be evaluated in terms of determinant formulas [47] Z(W ) = (±) det 2 cosh W 2 , {W 1 , W 2 } = (±) det(1 + e W 1 e W 2 ), (C. 4) where the ± in parentheses symbolise the eventual sign ambiguity. The square root (and hence the sign ambiguity) indicates that the pairs of eigenvalues of the skewsymmetric matrices must be taken into account with halved degeneracy [47]. Note that, for Gaussian states commuting with the particle number operator (i.e., when the exponent can be written with a Hermitian matrix in terms of the fermion operators instead of Majoranas), similar trace formulas apply without square roots and sign ambiguity [57]. In Section 5 we need the traces of operators O 3 + and O 2 + O − , respectively. Applying Eqs. (C.3) and (C.4), using hyperbolic identities for multiple arguments, one observes that the formulas can be expressed solely in terms of Γ ± with the result The sign ambiguity can be fixed by comparing to exact calculations of the traces. We find that the negative sign in the first trace of Eq. (C.5) is needed only in case Γ + contains a quadruplet of purely imaginary eigenvalues. For the Ising chain, this can happen only in the symmetry-broken phase, h < 1.
The numerics for small chains shows that, gradually decreasing the value of h, the appearance of this quadruplet exactly coincides with the vanishing of the first determinant. Interestingly, the second determinant in Eq. (C.5) always remains positive and thus no sign ambiguity appears there.