New Journal of Physics

The open access journal atthe forefront of physics

Deutsche Physikalische Gesellschaft DPG

IOP Institute of Physics

Published in partnership with: Deutsche Physikalische Gesellschaft and the Institute of Physics

OPEN ACCESS

RECEIVED

19 February 2015

REVISED

22 April2015

ACCEPTED FOR PUBLICATION

27 April2015

PUBLISHED

27 May 2015

Content from this work may be used under the terms ofthe Creative Commons Attribution 3.0 licence.

Any further distribution of this work must maintain attribution to the author(s) andthetitleof the work, journal citation andDOI.

On the partial transpose of fermionic Gaussian states

Viktor Eisler1 and Zoltan Zimboras2

1 MTA-ELTE Theoretical Physics Research Group, Eotvos Lorand University, Pazmany Peter setany 1/a, H-1117 Budapest, Hungary

2 Department of Computer Science, University College London, Gower Street, WC1E6BT London, UK

E-mail: zimboras@gmail.com

Keywords: quantum mechanics, Gaussian states, entanglement, fermion systems, negativity

Abstract

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 ofthe 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.

1. 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 ofthese systems are given bybosonic/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 well-behaved 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-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

© 2015 IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

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-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 section 3, focusing on a particular choice of basis which leads directly to our main result. Section 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 section 5. The paper concludes in section 6 with a short discussion ofthe results and their possible extensions. Various details ofanalytical calculations are included in three appendices.

2. Model and definitions

We consider quantum systems associated to free-fermion Hamiltonians

m, n=1

Amncm cn + ^ Bmncm cn

where the matrices A and B are Hermitian and antisymmetric, respectively. The fermionic creation/annihilation operators, cm and cm, satisfy the canonical anticommutation relations {cm, cn} = Smn. For our purposes, it will sometimes be more convenient to work with Majorana fermions defined as

a2m-1 cm + cm, a2m i( cm cm), (2)

satisfying the relations {ak, ai\ = 28ki. In terms ofMajorana operators, the Hamiltonian of equation (1) with

real A and B can be rewritten as H = Tm^man, where T2m,2n-1 = -?2n- 1,2m = 4(Amn + Bmn) and

T2m,2n = T2m-1,2n-1 = 0. The product of all Majorana operators define the parity operator, P = iN n^ an, 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-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-39]. A state p is Gaussian if it can be written as

2 Wkiaiax!4 k,l

where Wis a purely imaginary antisymmetric matrix (with the possibility of | Wki | ^ allowed) and Z is the normalization factor. A Gaussian state can also be characterized uniquely by its covariance matrix r = ([ ak, ai ])/2 via

tanh W = r. (4)

Using the covariance matrix, one can express the expectation value of any Majorana monomial through the Wick expansion

Tr(p aman2...an2e) = 2 sgn (n) H Fn

nn (2k—1), nn (2k)>

where all the indices are different, the sum runs over all pairings n, and sgn (n) denotes the sign of n 3. Let us note, that similarly to Gaussian states, one can introduce general Gaussian operators which are also defined through equation (3), however, without requiring that the spectrum of Wis real. The Wick expansion, i.e. equation (5), holds for these operators, as well.

We will also study spin chain models that are related to free-fermion Hamiltonians through the Jordan-Wignertransformation [40]

a 2m-1

H ak ^

where o^ (with a = x, y, z) denote the Pauli matrices acting on site m. The prototypical example is the transverse field Ising (TI) chain

hti = -_ ^ (xax +1 + hot),

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 equation ( 1 ) with

_ ( °m,n— 1 +

m, n— 1 + Om, n+1 )

3mn — 2 ( Om,n— 1 °m,n + 1j

Such a quadratic Hamiltonian can be diagonalized by a canonical transformation

and brought into the standard form

2 2[(m) a2m—1 — 'wk (m) a2m],

m— 1

H = 2Ak nl nk +

const.

k— 1

The spectrum Ak and the vectors yrk in equation (9) follow from the eigenvalue equations

(A - B)(A + B) $k = Ak fa,

(A + B)(A - B)Vk = Ak¥k.

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

(8) (9)

(11) (12)

with matrix elements given by

n21 n22

0 —ie

Smn = (m) (n)tanh

In our studies we will also be interested in the periodic TI chain, given by a Hamiltonian as in equation (7) with the sum running up to L and boundary condition o*+1 = o*. 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, (m) ~ elpkm and yk (m) ~ el6k elpk™, with the Bogoliubov angles and eigenvalues given by

h — eipk

Ak = ^ 1 + h2 — 2h cos pk ,

eiök =

Apairing n overaset {1,2,___, 2^} is a permutation of the 21 elements which satisfies n (2k — 1) < n (2k) and n (2k — 1) < n (2k + 1)

forany 1 ^ k ^ I. 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 (n) = (—1)M.

B A, B A2 B B A, A2 B

"(a)' * * * * * * * *

Figure 1. Two possible partitioning of a spin/fermionic chain into subsystems A1, A2, and B, as described in the text.

and the allowed values of the momenta are pk = (2k — 1) njL with k = —L/2 + 1,..., L/2. One can also work directly in the thermodynamic limit, L ^ x>, where the momentabecome continuous and the sum in equation (14) for the matrix elements gmn is replaced with an integral.

3. 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-43] and fermionic models [34,35,44,45] 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,46]. 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.

3.1. Definition of the partial transpose

Consider a general tripartition of a chain of qubits into disjoint sets A1, A2 and B, e.g. as in figure 1. Let p denote the density matrix of the whole composite system. Defining A = A1 U A2, the reduced density matrix (RDM) of subsystem A is given by pA = TrB p. The partial transpose of the RDM, pJ2, with respect to the subsystem A2 is defined by its matrix elements as

(ej pT2 |ek" e®) = (e/V | pA ^j, (16)

where {)} and {| ej2))} denote complete bases on the Hilbert spaces m1 and h2 pertaining to the subsets A1

and A2. The definition of pJ is basis dependent. However, one can easily characterize the set of transpositions on the operators acting on h2 as those non-degenerate linear transformations that satisfy

r (M1M2) = r (M2)r (M1), (17)

for any two operators M1 and M2 acting on h2. Since any two partial transpositions can be connected by a unitary conjugation, the eigenvalues of pJ2 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 ofa system with N fermionic modes, e.g. as in figure1.Let {mb m2, ..., m2k} and {n1, n2, ..., n2l} denote the indices of the Majorana operators belonging to the subsystems A1 and A2, respectively. Let us introduce the notation a0 = ! and a0 = ax. A general fermionic state on A = A1 U A2 can be written as

pa = 2wk,T < . a£kani . aZ, (18)

where the variables k = (k1, ... K2k) and t = (t1, ... , t21 ) in the summation run over all bit-strings of length 2k and 21, respectively. Note that since physical fermionic states must commute with the parity operator, as

discussed in section 2, one has wK,T =0when ^ ki + 2 1 T is odd.

The partial transpose of pA is simply a transformation that leaves the A1 Majorana modes invariant and acts as a transposition on the operators built up from modes of A2, i.e.

pA = 2Wk,T am\ . aKm2kr ( a£ ... a), (19)

where r satisfies equation (17). Since also in the fermionic case all the transpositions are connected by a unitary conjugation, the eigenvalues of pT will be independent ofwhich r we choose. It will be useful to consider the

following particular transposition which is defined by

' |i| mod 4 G {0, 3},

n (<...< ) = (-D/ (t) , where f (t) = j 0 f

t| mod 4 G {1, 2},

where |t | = 2-i T'. In other words, a Majorana monomial is mapped by r toitselfifitisoflength 4 n 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 [47]. Although the definition of entanglement in fermionic systems is somewhat different from the case of spin systems, it has been proven that pJ2 can only have negative eigenvalues if pA 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 A1 U A2 in a spin chain model and its fermionic counterpart are usually not equivalent (not isospectral), unless A1 and A2 are adjacent intervals, as depicted in figure 1 (b) [48,49]. 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.

3.2. The Gaussian case

Consider a Gaussian state pA on the system A = A1 U A2 with a covariance matrix

rA=(£! £2) (21)

where r11 and r22 denote the reduced covariance matrices of subsystems A1 and A2, respectively; while ru and r21 contain the expectation values of mixed quadratic terms. Let P2 be the parity operator on subsystem A2, and define the operators p+ = 1(pA + P2 pAP2) and p— = 1(pA — P2 pA P2).Bydefinition, we have that pA = p+ + p_. Using the notation of section 3.1, p+ and p_ can be expanded as

n = ^ w a K1 a K2k a T1 a T2l

F+ ~ ¿J K,T Wm1 ■■ ■ wmlkwn1 ■■ • "n2l >

It I even

p— = 2 wK,T < .•• . an%, (22)

|t| odd

where the coefficients wK,T can be obtained from rA using the Wick rule, equation (5). By linearity of the partial transpose, pj2 = p+2 + p-2 follows, and pj2 canbe obtained using equation (20):

p+T = 2 (-1)1

T |/2 w a K1 a K2k a T a

wK,l ami- - - am2k ani " - an2i,

K,T |t| even

2 (-1)(UI-1)2WKT ami...am*taZ...a%. (23)

ItI odd

Let us introduce the generalized Gaussian operators 0+ and O—, with covariance matrices

r+ = (ir21 — r22} r — = (— ir21 — r22} (24)

and consider the Majorana monomial expansion ofthese operators

O = 2n ± a K1 a K2k a T1 a T2l (25)

w± _ ¿J UK,T wm1---wmikwn1--- n2l'

Since 0+ and 0— are Gaussian operators, one can again obtain o±T from r± using equation (5). Connecting the Wick-expansion form of wKT with that of n±T, using the relation between rA and r±, one can deduce that

i±i( —1)(UI—1)/2wKT when |t| odd, o±t = i , .„ (26)

I ( —1pl/2wK,T when |t| even.

Comparing this with equation (23) it immediately follows that p+2 = 1(0+ + 0—) and p—2 = "^(0— — 0+). Thus, we obtain the decomposition

r 1 - i 1 + i

pT = — 0+ + — 0_, (27)

which is, from a conceptual point of view, the main result of the paper.

4. 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 0+ and 0_ do not commute in general, equation (27) can not be rewritten for the eigenvalues, and thus one does not have a direct access to the full spectrum of pj2. 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 ofentanglement.

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 equation (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, equation (7), with a small number of spins.

4.1. Bipartite pure states

We first consider the simplest case with B = 0 and a pure state on A = A1 U A2 given by p = |. Since for

any pure Gaussian fermionic state r2 = I is satisfied, it follows that [/+, r_] = 0, which implies that 0+ and 0_ commute as well. Furthermore, since their eigenvalues are connected by a complex conjugation, the spectrum of pT2 is simply given as the sum of the real and imaginary parts of the 0+ eigenvalues. Now, since 0+ is also Gaussian, its spectrum is uniquely determined through the eigenvalues of r+. To obtain them, we first assume without loss of generality that | A11 ^ | A21. Furthermore, for notational simplicity we also assume that | A | is even, however, the results for the odd case follow trivially. Let denote the eigenvalues of the reduced covariance matrix r11 with k = 1, ..., |A1| and Mk ^ 0. As shown in appendix B, the eigenvalues ±v± of r+ can then be given as

v± = J Mk ± _ m' when k = ^ IM (28)

k | 1 when k = |A1| + 1, ..., |A|/2.

The canonical diagonalized form of the Gaussian operator 0+ reads

A/2 T + inckb ck b ck

0+ = n n T + kl2k_1bk , (29)

k=1 ck=±

where b± are Majorana operators obtained from aj via the orthogonal transformation which diagonalizes r+. The eigenvalues of 0+ can be obtained according to the following rules. First, we shall consider the various combinations of the conjugate eigenvalue pairs for each k = 1, ..., |A1| as

1 + qk V+ 1 + o'v-= "!--!" ' (30)

-,OkOk

with ck = ± and c' = ±. Using equation (28) this yields

++ 1 + Mk __ 1 _ M +_ _+ in-2

mk =—2—, =—2—, m+ = _®k = 2V1 _ Mk , (31)

where we used the property v+v_ = 1. The nonzero eigenvalues of 0+ can then be written down as

|Al1 1 + okvk n ok ^^¡T^Mk

= n = n ^F1 n -^T^, (32)

Ok= ot Ok=-Ok

where o and O are the signature arrays corresponding to the eigenvalue. Note that the additional v± = 1 in equation (28) lead to further eigenvalues of 0+ that are all equal to zero.

The products in equation (32) are either real or purely imaginary and the eigenvalues of pT2 thus follow as Re Qss> or Im Qss', respectively. It is instructive, however, to derive the same spectrum using the Schmidt decomposition of |k). Dividing the Hilbert space h = h1 ® h2 into two parts4, one has

№> = k) №)> (33)

with the RDM eigenvalues X{ of subsystem A1. Clearly, the state is supported on a smaller Hilbert space H ® H and i ofp,

and is invariant under the action of a flip operation defined by | k/ )| kj) ^ I k j )l4't2). The partial transpose

PT2 = (I KK)T2 = |K) ® j K|, (34)

commutes also with the flip operator. Furthermore, it is easy to check that the eigenvalues and vectors are

pt2 K)K) = * i K) pT2 |#) = ±^ |k±), (35)

where we introduced the notation

\*±) = 7r(( j ± K ;) k» ^ j. (36)

Note that all the positive (negative) eigenvalues correspond to even (odd) eigenvectors with respect to the flip operation. Moreover, since p is Gaussian, one can immediately write down the products of eigenvalues as

j— n S, (37)

dt= a' Ck=- a'

where the signature a has again components ak = ±. Comparing equation (37) to ( 32), one indeed recognizes Qsa> up to the factors of i.

Owing to the simple product structure of the eigenvalues in equation (37) and, in particular, to the fact that all the negative eigenvalues are located in the odd subspace, one has

Tr \pT\ = 1 - 2Tro pT = H ( + V1 - M2 ), (38)

where we introduced the notation Tre/o for traces taken over the even/odd subspace. Thus, the pure-state logarithmic negativity is given by the simple formula

£ = ln Tr |pT2| = 2 ln(l + 71 - )• (39)

Considering also the expression of the Rényi entropy for fermionic Gaussian states

1 +1L j" + | 1 - ^

S = ? - [( ^ )" + ( ^ j"

the well-known equality £ = Si/2 for pure states can be confirmed directly.

4.2. 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 w.r.t. 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 |Ai | = |A21 = N/2 and states that are reflection symmetric. To analyse such a situation, we shall consider Gibbs states of the open TI chain, equation (7), 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.

Considering a covariance matrix of the form (13), the spectrum of r+ is invariant with respect to a sign change and complex conjugation, hence the eigenvalues can be collected into two families of quadruplets

{, z* -zk, -zk}, k e (I) {{, -ivk, -iuk, ivk}, k e (II), (41)

When considering a tensor product structure, we always refer to a partition of the spin chain system constructed through the Jordan-Wigner transformation.

where in family (I) we choose Re Zk > 0 and Im Zk > 0,whereas Uk > Vk infamily (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 Uk # Vk. 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 equation (7).

Analogously to the pure case in equation (30), we first assign the factors

fflf Ok = ^

^ Ok =

1 + |zk I2 + Ok2 Re Zk ,

---> Ok = °k

1 - |zk |2 + Ok2i Im Zk

, k e (I)

1 + UkVk + Oki(u - v) 4

1 - UkVk + Oki(u + v)

Ok = -Ok

, Ok = Ok Ok = -Ok

, k e (II) (42)

within each quadruplet, and the eigenvalues QO o of 0+ are again given in the factorized form of equation (32). Although the spectrum of the operator 0- is identical to that of 0+, they do not commute in general and thus one has no direct access to the eigenvalues of p12. Nevertheless, the information about the even/odd parity of the eigenvectors is retained. In fact, the reflection operator R, which defines the even/odd subspaces in our case, acts on the spin operators as Ro^Rt = o^-n (with a = x, y, z), implying the action RcjRt = PcN-n on the creation operators, where P is the parity operator. Using this, it follows that the sign factor associated to an eigenvector of 0+ reads

n-.—11, if Ok = ok,

Sok o' + Im M Sok o^ Sok o'=< . ., , (43)

k k I Okh if Ok = -Ok,

which can also be verified by considering the pure state limit. The parity of the 0- eigenvectors are simply obtained through the factors s*O'.

With the knowledge of SO ¿, we are now able to carry out signed traces of the form

Tre0+ - Tro0+ = £soOQoo-- (44)

Note that the terms in the sum of equation (44) completely factorize in the quadruplet index k. Furthermore, using the fact that Tre/o0_ = (Tre/o0+)* and thus Tre/opr2 = ReTre/o0+ + ImTre/o0+, a simple calculation leads to

T r2 T -pr 1 + |zk|2 + 2Imzk Y-r 1 + ukvk + uk + vk

1 - 2Tro p12 = Tre p12 - Tro p12 = n -2- n -2-' (45)

ke(I) 2 ke(II) 2

Finally, we define the quantity

So = ln max( 1 - 2 Trop12, 1), (46)

which clearly gives a lower bound for the logarithmic negativity, 8 ^ 8^, with strict equality if all the negative eigenvalues reside in the odd sector and their number is equal to the dimension ofthat 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 [50,51].

To test the bound 8o against the exact value of the logarithmic negativity, we considered small TI chains with N ^ 10 and obtained 8 via exact diagonalization of p12. This is shown on figure 2, as a function of the temperature (left) as well as of the magnetic field (right). We find that, for low enough temperatures, 8o indeed exactly coincides with 8. 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 ofthese processes increasing the difference 8 - 8o. Nevertheless, for the small system sizes considered, it appears that 8o gives a very good approximation up to temperatures T « 1/N, above which it starts to deviate significantly.

4.3. 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 figure 1 (b). In this case, the reduced state pA after tracing out the sites of B is a mixed Gaussian state, associated to the reduced covariance matrix rA, with indices running over sites in A [52].

Figure 2. Logarithmic negativity 8 (symbols) versus 8o (lines) between two halves of a small Ising chain with Hamiltonian (7)ina Gibbs state. Left: as a function of the inverse temperature f with h = 1 and various number of spins N. The curves 8 and 8o start to deviate around /) » N. Right: as a function of h for various f and N = 8. Minor deviations between 8 and 8o are only visible for f = 10.

The logarithmic negativity for the tripartite case can be obtained with CFT methods [23]. For two intervals of the same size I embedded in a system of length L with periodic boundary conditions, the calculation yields [24]

£ (I, L) = - ln 4

n . (I )

sin (T )

+ const.,

with the central charge c and a non-universal constant. However, subtracting the value at I = LI4 one obtains a universal scaling function

e (z) = £ (I, L) - £ (LI4, L) = 1 ln [tan(nz)],

where z = IL 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 equation (46), for the geometry at hand. In fact, since reflection symmetry is fulfilled, all the arguments of the previous section, leading to equation (45), apply and 8o is given by the same formula in terms of the eigenvalues of r+. However, there is an important difference compared to the bipartite thermal case, which is apparent from the numerical investigation ofsmall systems. Namely, the number ofnegative eigenvalues of pAT2 is always less then the dimension ofthe odd subspace and, moreover, some of the corresponding eigenvectors are even. Thus, by tracing out the sites of B, the partial transpose pj2 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 8o is tight.

In spite of the above findings, 8o shows a very interesting behaviour, which is demonstrated on figure 3. First of all, it shows a clear signature of the phase transition at h = 1, which can be seen when plotting 8o (LI 4, L) against h on the left of figure 3. Furthermore, defining the quantity eo analogously to equation (48), one finds an excellent data collapse when plotted against the variable z = f/L, see right of figure 3. The scaling function is found to be given by

eo(z) = £o(l, L) - £0(L/4, L) = — ln

tan(nz) 2 cos2(nz)

Although the functional form of eo (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 e (z) in equation (48), however, the argument is modified as well. We also performed a calculation directly in the thermodynamic limit L ^ rc>,andfound 8o(l) = 1/16 ln I + 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 Bmn = 0 and Amn = ^(¿im, n-1 + <5m, n+1). In complete analogy with the result for the

bipartite entanglement [53], one finds 8XX (21) = 28TV) for h =1, and thus a doubled prefactor 1/8 with respect to the critical TI chain. Therefore, even though 8o does not approximate 8 well, it shows exactly the same universal behaviour, suggesting it as an entanglement indicator which is extremely simple to calculate.

Figure 3. The scaling of £o (I, L) for two adjacent intervals of equal length I in the ground state of a periodic Ising chain with L sites and magnetic field h. Left: £o (LI 4, L) as a function of the magnetic field h for various L. The logarithmic divergence around h =1is clearly seen. Right: £o(l, L) — £o(LI 4, L) at the critical point h = 1 against the scaling variable z. The solid line shows the scaling function in equation (49).

5. Trace formulas

Our result for the partial transpose, equation (27), can be further tested by looking at traces of integer powers of the partial transpose pj2 which can also be carried out within CFT. Since the identity Tr(pJ2 )2 = Tr pA 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]

Rs(l, L) = lnTr(pT)3 = — 9 ln L- sin2 (() sin (( Similarly to equation (48), a universal scaling function can be defined as [24]

r3(z) = Rj(l, L) — R3(L!4, L) = — - ln [2 sin2(nz)sin(2nz)],

+ const.

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 I, embedded in an infinite chain which is thermalized at inverse temperature Applying a simple conformal transformation, the corresponding CFT formula follows as

R3(l, P)

P3sinh2 ( nL J sinh (

+ const.

We now show how the above traces can be calculated with the covariance matrix formalism. Expanding the

third power of p T in equation (27) and taking the trace one arrives at

Tr(pT)3 = -2 Tr( 0+ ) + |Tr(0+ 0-),

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 [48,54]. The main steps of this calculation are summarized in appendix C. In turn, one finds

Tr(PÎ)3

1 + r+ + 2r+r - j

where the sign of the first term depends on the spectrum of /+, see equation (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 [48]. 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 r± and evaluating the determinants. This is shown in figure 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.

Figure 4. Left: R3 (1, L) — R3 (LI 4, L) as function of z = l/L for two adjacent intervals of length I in the ground state of the critical TI chain of size L. The solid line shows the CFT formula (51). Right: R3 (1, f) for adjacent intervals of size I in a thermal state of the infinite chain (h = 1), with inverse temperature/). The inset shows the rescaled data compared to the CFT prediction (52) on a horizontal log-scale.

6. 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 r±. 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 pT 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 £o with earlier studies [55-57].

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 [48]. 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 ofthe traces for such a case is still realizable, but presumablymore tedious.

It would also be interesting to see whether the lower bound £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 equation (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 £o and £ 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.

Finally, we point out that some specific classes of mixed Gaussian states exist which allow for an exact calculation of the logarithmic negativity using the methods introduced here. These are the states for which the relation [/+, r—] = 0 is satisfied, an example being the isotropic Gaussian states [44], for which the covariance matrix satisfies r2 = X2i with some 0 ^ X < 1. The situation is then analogous to the pure-state case and the corresponding calculation can be generalized, which we leave for future studies.

Acknowledgments

We would like to thank Leonardo Banchi and Fernando Brandao for useful discussions, and an anonymous referee for pointing out [44]. 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 ofthe Gaussian RDM is given by

1 + iVkb 2k —1 b 2k

where Vk G [0, 1] and the Majorana operators bj are related to aj through an orthogonal transformation. For simplicity, we will consider only covariance matrices of the form of equation (13), and assume also reflection symmetry. These states are parametrized by a single angle в beside the covariance matrix eigenvalues Vk.

Using the Jordan-Wigner representation of aj and working in the usual spin basis, the most general form of the partial transpose is

1 V1V2

p12 = - +

-öl Ö2 +

Vi + V2

cos 2Q

0 sin 2в sin 2в 0 0 -cos 2в

V1 - V2 0 0

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 + V1V2 ± V(Vi + V2)2 cos2 29 + (Vi — V2)2 ]/4,

A3,4 = [ 1 — v1v2 ± ( v1 + v2)sin 29]/4.

Using the parity operator P = <72z, one can also construct p±2 = (pT ± P2 pTl P)j2 with matrix elements

cos 29 0

0 0 0 0

, 0 —cos 29,

T2 1 V1V2 z z ^ V1 + V2

P+-2 = 7 + _T"ö1 ö2 + ---

+ 4 4 4

T2 V1 + V2

P-2 = —T~

0 sin 2в sin 2в 0 0 0

V1 - V2 0 0

The eigenvalues of the operator p+2 + ip—2 then read

Q1,2 = [ 1 + V1V2 (V1 + V2)2 cos2 29 — (V1 — V2)2 ]/4,

ß3,4 = [ 1 - V1V2 ± i( v1 + v2)sin 2в]/4.

Note that we have Д3)4 = Re ^3 4 + Im ^3 4 which simply follows from the fact that p+2 and p—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 0+ with covariance matrix Г+ has indeed eigenvalues given by equation (A.6). The covariance matrix Г for the Gaussian state (A.1) and the associated Г+ have the form

Г = i

0 c 0 \ s- ' 0 c 0 \ is-

- c 0 s+ 0 - c 0 is+ 0

0 0 , Г+ = i

-s+ c 0 is+ 0 - c

-s- 0 -c 0 -is- 0 c 0

with the shorthand notation

V1 + V2 c = -cos 2в,

2 2 The four eigenvalues ± v± of Г+ can be computed with

V1 + V^ V1 - V2

s± =-sin 2в ±-.

S+ + s-

(A. 9)

If the operator O+ is Gaussian, its eigenvalues must have the form

^++ = i±v+, fl„ 1 - v+1

2 2 2 2

1 + V+ 1 - V- 1 - V+ 1 + V-

Q+- =-+-, Q-+ =-+-. (A.10)

+ 2 2+ 22

Substituting (A.8) and (A.9) into (A.10), we indeed recover the values in (A.6).

Finally, let us shortly discuss the non-Gaussian character of pTl. comparing the expectation values Tr(pr2 aman),where m, n = 1, ..., 4,with Tr(pr2 a1a2 a3 a4), one observes that the Wick expansion, equation (5), does not hold, unless c2 = v1v2. 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.

Appendix B. Eigenvalues of r± for pure states

Here we consider the eigenvalue problem of the modified covariance matrices r±, that are associated to a pure-state covariance matrix r of the form (13). The RDMsfor subsystems A1 and A2 are determined via the reduced covariance matrices r11 and r22. Since they split into two submatrices, one could equivalently consider the squared eigenvalue problem of the matrices

T-i T-i lv (m) m G Aa,

GZ = 2gmlgnl = YMnv;(m)<(n), (m) = j \ (B.1)

I 0 m G Aa,

lGa p,q L 1-

with nonzero matrix elements only within m, n G Aa ,a = 1, 2. The overlap matrices Ma have matrix elements

f (l) l G Aa,

mm= xk(l)<(l), <(l) = \ \ lctA (b.2)

l U l £ Aa.

Notethatboth G11and G22 have the same eigenvalues < 1 with k = 1,..., min(|A1|, |A21 ), whereas ^k = 1 for the remaining eigenvalues of the larger matrix. We also introduce the block-diagonal matrix

G = iG11 (B.3)

i 0 G22 J

with all the nontrivial eigenvalues being doubly degenerate.

The matrix elements of the covariance matrices r± in equation (24) are determined through

g±n = 2v± (m) 4>± (n), (B.4)

with the vectors

(l) = ^ (l) ± i^q2 (l), (l) = V\ (l) ± iVq2 (l). (B.5)

Thus the spectrum of r± follows from the eigenvalues of the squared matrix

G±n = 2 (Mpq - MPq )± (m) Vq± (n). (B.6)

Inserting (B.5) and using the completeness property Mlpq + Mpq = Spq, the matrices G± have the block form

G±= i2G11 "r 1 ±f ), (B.7) i ±2iFT 2G22 - 1 )

Fmn =2Mjq Vp (m) Vq(n). (B.8)

It is easy to check that F satisfies

FFT = G11( 1 - G11), FTF = G22( 1 - G22), (B.9) and thus the following matrix identity holds

( G±- 2 G + 1)2 =-4G (1 - G). (B.10)

Rewriting in terms of the eigenvalues (vj±)2 and pk of G± and G, respectively, one arrives at

( v±)2 = 2^k2 -1 ± 1 - m2 = ( ± W1 - M2 )2.

(B.11)

Appendix C. Determinant formulas

Let us consider the Gaussian operators 0± corresponding to the generalized covariance matrices r± = tanh (W±/2). With a denoting the vector of Majorana operators, we introduce

and thus 0± = 0 (w±)Z (W±) Our aim is to calculate various traces of the form Tr( 0+n 0") with some integers m and n. Following the lines of [48], we first introduce the notation

{W W2} = Tr [0 (W) (W2)]. We also note the simple fact that 0m (W±) = 0 (mW±). Hence, the traces we consider can be

written in the form

{mW+, "W-}

( o+mo-) = —--

They can be evaluated in terms of determinant formulas [48]

Z(W) = (±) Jdet(2 cosh — j, {W1, W2} = (±)^det( 1 + eW1eW2),

Z (W) = (±)

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 skew-symmetric matrices must be taken into account with halved degeneracy [48]. 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 [58].

In section 5 we need the traces of operators 0+ and 0+ 0—, respectively. Applying equations (C.3) and (C.4), using hyperbolic identities for multiple arguments, one observes that the formulas can be expressed solely in terms of r± 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 equation (C.5) is needed only in case r+ contains a quadruplet of purely imaginary eigenvalues. For the Ising chain, this can happen onlyin 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 equation (C.5) always remains positive and thus no sign ambiguity appears there.

References

[1] Amico L, Fazio R, Osterloh A and Vedral V 2008 Rev. Mod. Phys. 80 517

[2] Calabrese P, Cardy J and Doyon B 2009 J. Phys. A: Math. Theor. 42 500301

[3] Eisert J, Cramer M and Plenio M B 2010 Rev. Mod. Phys. 82 277

[4] Vidal G, Latorre J I, Rico E and Kitaev A 2003 Phys. Rev. Lett. 90 227902

[5] Calabrese P and Cardy J 2009 J. Phys. A: Math. Theor. 42 504005

[6] Peschel I and Eisler V 2009 J. Phys. A: Math. Theor. 42 504003

[7] Plenio M and Virmani S 2007 Quantum Inf. Comput. 7 1

[8] Peres A1996 Phys. Rev. Lett. 77 1413

[9] Horodecki M, Horodecki P and Horodecki R1996 Phys. Lett. A 223 1

[10] Vidal G and Werner R F 2002 Phys. Rev. A 65 032314

[11] Plenio M B 2005 Phys. Rev. Lett. 95 090503

[12] Simon R 2000 Phys. Rev. Lett. 84 2726

[13] Audenaert K, Eisert J, Plenio M B and Werner R F 2002 Phys. Rev. A 66 042327

[14] Ferraro A, Cavalcanti D, Garcia-Saez A and Acin A 2008 Phys. Rev. Lett. 100 080502

[15] Cavalcanti D, Ferraro A, Garcia-Saez A and Acin A 2008 Phys. Rev. A 78 012335

[16] Anders J and Winter A 2008 Quantum Inf. Comput. 8 0245

[17] Anders J 2008 Phys. Rev. A 77 062102

[18] Marcovitch S, Retzker A, Plenio M B and Reznik B 2009 Phys. Rev. A 80 012325

[19] Wichterich H, Molina-Vilaplana J and Bose S 2009 Phys. Rev. A 80 010304(R)

[20] Bayat A, Sodano P and Bose S 2010 Phys. Rev. B 81 064429

[21] Wichterich H,Vidal J and Bose S 2010 Phys.Rev. A 81 032311

[22] Santos R A, Korepin V and Bose S 2011 Phys. Rev. A 84 062307

[23] Calabrese P, Cardy J and Tonni E 2012 Phys. Rev. Lett. 109 130502

[24] Calabrese P, Cardy J and Tonni E 2013 J. Stat. Mech. P02008

[25] Nobili C D, Coser A and Tonni E 2015 arXiv: 1501.04311

[26] Calabrese P,CardyJ and Tonni E 2015 J. Phys. A: Math. Theor. 48 015006

[27] Eisler V and Zimboras Z 2014 New J. Phys. 16123020

[28] Coser A, Tonni E and Calabrese P 2014 J. Stat. Mech. P12017

[29] Hoogeveen M and Doyon B 2014 arXiv:1412.7568

[30] Wen X, Chang P and Ryu S 2015 arXiv: 1501.00568

[31] Calabrese P, Tagliacozzo L and Tonni E 2013 J.Stat. Mech. P05002

[32] Alba V 2013 J. Stat. Mech. P05013

[33] Chung C, EAlba V, Bonnes L, Chen P and Lauchli A M 2014 Phys. Rev. B 90 064401

[34] Banuls M C, Cirac JI and Wolf M M 2007 Phys. Rev. A 76 022311

[35] Banuls M C, Cirac JI and Wolf M M 2009 J. Phys. Conf. Ser. 171 012032

[36] Zimboras Z,ZeierR,KeylMandSchulte-HerbruggenT2014EPJQuantum Technol. 1 11

[37] Bravyi S 2005 Quantum Inf. Comput. 5 216

[38] BravyiSandKonigR2012 Quantum Inf. Comput. 12 925

[39] deMeloF,CwiklinskiPandTerhalBM2013NewJ.Phys. 15013015

[40] Lieb E, Schultz T and Mattis D 1961 Ann. Phys., NY 16 407

[41] Werner R F and Wolf M M 2001 Phys. Rev. Lett. 86 3658

[42] Giedke G, Kraus B, Lewenstein M and Cirac J12001 Phys. Rev. Lett. 87167904

[43] Weedbrook C, Pirandola S, Garcia-Patron R, Cerf N J, Ralph T C, Shapiro J H and Lloyd S 2012 Rev. Mod. Phys. 84 621

[44] Botero A and Reznik B 2004 Phys. Lett. A 331 39

[45] Benatti F, Floreanini RandMarzolino U 2014 Phys. Rev. A89 032326

[46] Serafini A, Adesso G and Illuminati F 2005 Phys. Rev. A 71 032349

[47] Sarosi G and LevayP 2014 J. Phys. A:Math. Theor. 47115304

[48] Fagotti M and Calabrese P 2010 J. Stat. Mech. P04016

[49] Igloi F and Peschel 12010 Europhys. Lett. 89 40001

[50] Johnston N2013 Phys. Rev. A 87 064302

[51] RanaS2013 Phys.Rev. A 87 054301

[52] Latorre JI and Riera A 2009 J. Phys. A: Math. Theor. 42 504002

[53] Igloi F and Juhasz R 2008 Europhys. Lett. 81 57003

[54] Banchi L, Giorda P and Zanardi P 2014 Phys. Rev. E 89 022102

[55] TothG2005 Phys.Rev. A 71010301

[56] Wu L A, Bandyopadhyay S, Sarandy M S and Lidar D A 2005 Phys. Rev. A 72 032309

[57] Guhne O and Toth G 2009 Phys. Rep. 474 1

[58] Klich 12003 Quantum Noisein MesoscopicPhysics (NATO Science Series IIvol 97) (Dordrecht: Kluwer) pp 397-402

Copyright of New Journal of Physics is the property of IOP Publishing and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.