On the partial transpose of fermionic Gaussian states

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 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 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 of the results and their possible extensions. Various details of analytical calculations are included in three appendices.

We consider quantum systems associated to free-fermion Hamiltonians

$$\begin{eqnarray}&&H=\mathop{\sum }\limits_{m,n=1}^{N}\left[ {{A}_{mn}}c_{m}^{\dagger }{{c}_{n}}+\frac{1}{2}{{B}_{mn}}c_{m}^{\dagger }c_{n}^{\dagger }-\frac{1}{2}B_{mn}^{*}{{c}_{m}}{{c}_{n}} \right],\end{eqnarray} \tag{ 1 }$$

where the matrices A and B are Hermitian and antisymmetric, respectively. The fermionic creation/annihilation operators, $c_{m}^{\dagger }$ and c_m, satisfy the canonical anticommutation relations $\left\{ c_{m}^{\dagger },{{c}_{n}} \right\}={{\delta }_{mn}}$. For our purposes, it will sometimes be more convenient to work with Majorana fermions defined as

$$\begin{eqnarray}&&{{a}_{2m-1}}={{c}_{m}}+c_{m}^{\dagger },\qquad {{a}_{2m}}={\rm i}({{c}_{m}}-c_{m}^{\dagger }),\end{eqnarray} \tag{ 2 }$$

satisfying the relations $\left\{ {{a}_{k}},{{a}_{l}} \right\}=2{{\delta }_{kl}}$. In terms of Majorana operators, the Hamiltonian of equation (1) with real A and B can be rewritten as $H={\rm i}\sum _{m,n=1}^{2N}{{T}_{m,n}}{{a}_{m}}{{a}_{n}}$, where ${{T}_{2m,2n-1}}=-{{T}_{2n-1,2m}}=\frac{1}{4}({{A}_{mn}}+{{B}_{mn}})$ and ${{T}_{2m,2n}}={{T}_{2m-1,2n-1}}=0$. The product of all Majorana operators define the parity operator, $P={{{\rm i}}^{N}}\prod _{n=1}^{2N}{{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–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 ρ is Gaussian if it can be written as

$$\begin{eqnarray}&&\rho =\frac{1}{Z}{\rm exp} (\mathop{\sum }\limits_{k,l}\;{{W}_{kl}}{{a}_{k}}{{a}_{l}}/4),\end{eqnarray} \tag{ 3 }$$

where W is a purely imaginary antisymmetric matrix (with the possibility of $|{{W}_{kl}}|\to \infty$ allowed) and Z is the normalization factor. A Gaussian state can also be characterized uniquely by its covariance matrix ${{\Gamma }_{kl}}=\langle \left[ {{a}_{k}},{{a}_{l}} \right]\rangle /2$ via

$$\begin{eqnarray}&&{\rm tanh} \frac{W}{2}=\Gamma .\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&{\rm Tr}(\rho \;{{a}_{{{n}_{1}}}}{{a}_{{{n}_{2}}}}\ldots {{a}_{{{n}_{2\ell }}}})=\mathop{\sum }\limits_{\pi }{\rm sgn} \;(\pi )\mathop{\prod }\limits_{k=1}^{\ell }{{\Gamma }_{{{n}_{\pi (2k-1)}},{{n}_{\pi (2k)}}}},\end{eqnarray} \tag{ 5 }$$

where all the indices are different, the sum runs over all pairings π, and ${\rm sgn} \;(\pi )$ denotes the sign of π ³ . 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 W is 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–Wigner transformation [40]

$$\begin{eqnarray}&&{{a}_{2m-1}}=\mathop{\prod }\limits_{k=1}^{m-1}\sigma _{k}^{z}\;\sigma _{m}^{x},\quad \quad {{a}_{2m}}=\mathop{\prod }\limits_{k=1}^{m-1}\sigma _{k}^{z}\;\sigma _{m}^{y},\end{eqnarray} \tag{ 6 }$$

where $\sigma _{m}^{\alpha }$ (with $\alpha =x,y,z$) denote the Pauli matrices acting on site m. The prototypical example is the transverse field Ising (TI) chain

$$\begin{eqnarray}&&{{H}_{{\rm TI}}}=-\frac{1}{2}\mathop{\sum }\limits_{m}\left( \sigma _{m}^{x}\sigma _{m+1}^{x}+h\sigma _{m}^{z} \right),\end{eqnarray} \tag{ 7 }$$

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

$$\begin{eqnarray}&&{{A}_{mn}}=\frac{1}{2}({{\delta }_{m,n-1}}+{{\delta }_{m,n+1}})-h{{\delta }_{m,n}},\qquad {{B}_{mn}}=\frac{1}{2}({{\delta }_{m,n-1}}-{{\delta }_{m,n+1}}).\end{eqnarray} \tag{ 8 }$$

Such a quadratic Hamiltonian can be diagonalized by a canonical transformation

$$\begin{eqnarray}&&{{\eta }_{k}}=\mathop{\sum }\limits_{m=1}^{N}\frac{1}{2}\left[ {{\phi }_{k}}(m)\;{{a}_{2m-1}}-{\rm i}{{\psi }_{k}}(m)\;{{a}_{2m}} \right],\end{eqnarray} \tag{ 9 }$$

and brought into the standard form

$$\begin{eqnarray}&&H=\sum k=1N{{\Lambda }_{k}}\eta _{k}^{\dagger }{{\eta }_{k}}+{\rm const}.\end{eqnarray} \tag{ 10 }$$

The spectrum ${{\Lambda }_{k}}$ and the vectors ${{\phi }_{k}},{{\psi }_{k}}$ in equation (9) follow from the eigenvalue equations

$$\begin{eqnarray}&&(A-B)(A+B){{\phi }_{k}}=\Lambda _{k}^{2}{{\phi }_{k}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&(A+B)(A-B){{\psi }_{k}}=\Lambda _{k}^{2}{{\psi }_{k}}.\end{eqnarray} \tag{ 12 }$$

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

$$\begin{eqnarray}\Gamma =\left[ \begin{array}{ccccccccccccccc} {{\Pi }_{11}} & {{\Pi }_{12}} & \cdots & {{\Pi }_{1N}} \\ {{\Pi }_{21}} & {{\Pi }_{22}} & {} & \vdots \\ \vdots & {} & \ddots & \vdots \\ {{\Pi }_{N1}} & \cdots & \cdots & {{\Pi }_{NN}} \\ \end{array} \right],\qquad {{\Pi }_{mn}}=\left( \begin{array}{ccccccccccccccc} 0 & -{\rm i}{{g}_{nm}} \\ {\rm i}{{g}_{mn}} & 0 \\ \end{array} \right),\end{eqnarray} \tag{ 13 }$$

with matrix elements given by

$$\begin{eqnarray}&&{{g}_{mn}}=\mathop{\sum }\limits_{k}{{\psi }_{k}}(m){{\phi }_{k}}(n){\rm tanh} \frac{\beta {{\Lambda }_{k}}}{2}.\end{eqnarray} \tag{ 14 }$$

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 $\sigma _{L+1}^{x}=\sigma _{1}^{x}$. 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, ${{\phi }_{k}}(m)\sim {{{\rm e}}^{{\rm i}{{p}_{k}}m}}$ and ${{\psi }_{k}}(m)\sim {{{\rm e}}^{{\rm i}{{\theta }_{k}}}}{{{\rm e}}^{{\rm i}{{p}_{k}}m}}$, with the Bogoliubov angles and eigenvalues given by

$$\begin{eqnarray}&&{{{\rm e}}^{{\rm i}{{\theta }_{k}}}}=\frac{h-{{{\rm e}}^{{\rm i}{{p}_{k}}}}}{{{\Lambda }_{k}}},\qquad {{\Lambda }_{k}}=\sqrt{1+{{h}^{2}}-2h{\rm cos} {{p}_{k}}},\end{eqnarray} \tag{ 15 }$$

and the allowed values of the momenta are ${{p}_{k}}=(2k-1)\pi /L$ with $k=-L/2+1,...,L/2$. One can also work directly in the thermodynamic limit, $L\to \infty$, where the momenta become continuous and the sum in equation (14) for the matrix elements g_mn is replaced with an integral.

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 A₁, A₂ and B, e.g. as in figure 1. Let ρ denote the density matrix of the whole composite system. Defining $A={{A}_{1}}\cup {{A}_{2}}$, the reduced density matrix (RDM) of subsystem A is given by ${{\rho }_{A}}={\rm T}{{{\rm r}}_{B}}\;\rho$. The partial transpose of the RDM, $\rho _{A}^{{{T}_{2}}}$, with respect to the subsystem A₂ is defined by its matrix elements as

$$\begin{eqnarray}&&\langle e_{i}^{(1)}e_{j}^{(2)}|\rho _{A}^{{{T}_{2}}}|e_{k}^{(1)}e_{l}^{(2)}\rangle =\langle e_{i}^{(1)}e_{l}^{(2)}|{{\rho }_{A}}|e_{k}^{(1)}e_{j}^{(2)}\rangle ,\end{eqnarray} \tag{ 16 }$$

where $\{|e_{i}^{(1)}\rangle \}$ and $\{|e_{j}^{(2)}\rangle \}$ denote complete bases on the Hilbert spaces ${{\mathcal{H}}_{1}}$ and ${{\mathcal{H}}_{2}}$ pertaining to the subsets A₁ and A₂. The definition of $\rho _{A}^{{{T}_{2}}}$ is basis dependent. However, one can easily characterize the set of transpositions on the operators acting on ${{\mathcal{H}}_{2}}$ as those non-degenerate linear transformations that satisfy

$$\begin{eqnarray}&&\mathcal{R}({{M}_{1}}{{M}_{2}})=\mathcal{R}({{M}_{2}})\mathcal{R}({{M}_{1}}),\end{eqnarray} \tag{ 17 }$$

for any two operators M₁ and M₂ acting on ${{\mathcal{H}}_{2}}$. Since any two partial transpositions can be connected by a unitary conjugation, the eigenvalues of $\rho _{A}^{{{T}_{2}}}$ 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].

Zoom In Zoom Out Reset image size

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 figure 1 . Let $\{{{m}_{1}},{{m}_{2}},\ldots ,{{m}_{2k}}\}$ and $\{{{n}_{1}},{{n}_{2}},\ldots ,{{n}_{2\ell }}\}$ denote the indices of the Majorana operators belonging to the subsystems A₁ and A₂, respectively. Let us introduce the notation $a_{x}^{0}=$ and $a_{x}^{1}={{a}_{x}}$. A general fermionic state on $A={{A}_{1}}\cup {{A}_{2}}$ can be written as

$$\begin{eqnarray}&&{{\rho }_{A}}=\mathop{\sum }\limits_{\underline{\kappa },\underline{\tau }}{{w}_{\underline{\kappa },\underline{\tau }}}\;a_{{{m}_{1}}}^{{{\kappa }_{1}}}\;\ldots \;a_{{{m}_{2k}}}^{{{\kappa }_{2k}}}a_{{{n}_{1}}}^{{{\tau }_{1}}}\;\ldots \;a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}},\end{eqnarray} \tag{ 18 }$$

where the variables $\underline{\kappa }=({{\kappa }_{1}},\;\ldots \;{{\kappa }_{2k}})$ and $\underline{\tau }=({{\tau }_{1}},\;\ldots \;,{{\tau }_{2\ell }})$ in the summation run over all bit-strings of length $2k$ and $2\ell$, respectively. Note that since physical fermionic states must commute with the parity operator, as discussed in section 2, one has ${{w}_{\underline{\kappa },\underline{\tau }}}$= 0 when $\sum _{i=1}^{2k}{{\kappa }_{i}}+\sum _{j=1}^{2\ell }{{\tau }_{j}}$ is odd.

The partial transpose of ${{\rho }_{A}}$ is simply a transformation that leaves the A₁ Majorana modes invariant and acts as a transposition on the operators built up from modes of A₂, i.e.

$$\begin{eqnarray}&&\rho _{A}^{{{T}_{2}}}=\mathop{\sum }\limits_{\underline{\kappa },\underline{\tau }}{{w}_{\underline{\kappa },\underline{\tau }}}\;a_{{{m}_{1}}}^{{{\kappa }_{1}}}\;\ldots \;a_{{{m}_{2k}}}^{{{\kappa }_{2k}}}\mathcal{R}(a_{{{n}_{1}}}^{{{\tau }_{1}}}\;\ldots \;a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}}),\end{eqnarray} \tag{ 19 }$$

where $\mathcal{R}$ satisfies equation (17). Since also in the fermionic case all the transpositions are connected by a unitary conjugation, the eigenvalues of $\rho _{A}^{{{T}_{2}}}$ will be independent of which $\mathcal{R}$ we choose. It will be useful to consider the following particular transposition which is defined by

$$\begin{eqnarray}\mathcal{R}(a_{{{n}_{1}}}^{{{\tau }_{1}}}\ldots a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}})={{(-1)}^{f(\underline{\tau })}}a_{{{n}_{1}}}^{{{\tau }_{1}}}\ldots a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}},\ \ {\rm where}\ \ f(\underline{\tau })=\left\{ \begin{array}{ccccccccccccccc} 0 & {\rm if}|\underline{\tau }|\ {\rm mod} \ 4\in \{0,3\}, \\ 1 & {\rm if}|\underline{\tau }|\ {\rm mod} \ 4\in \{1,2\}, \\ \end{array} \right.\end{eqnarray} \tag{ 20 }$$

where $|\underline{\tau }|=\sum _{i=1}^{2\ell }{{\tau }_{i}}$. In other words, a Majorana monomial is mapped by $\mathcal{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 [47]. Although the definition of entanglement in fermionic systems is somewhat different from the case of spin systems, it has been proven that $\rho _{A}^{{{T}_{2}}}$ can only have negative eigenvalues if ${{\rho }_{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}}\cup {{A}_{2}}$ in a spin chain model and its fermionic counterpart are usually not equivalent (not isospectral), unless A₁ and A₂ 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 ${{\rho }_{A}}$ on the system $A={{A}_{1}}\cup {{A}_{2}}$ with a covariance matrix

$$\begin{eqnarray}{{\Gamma }_{A}}=\left( \begin{array}{ccccccccccccccc} {{\Gamma }^{11}} & {{\Gamma }^{12}} \\ {{\Gamma }^{21}} & {{\Gamma }^{22}} \\ \end{array} \right),\end{eqnarray} \tag{ 21 }$$

where ${{\Gamma }^{11}}$ and ${{\Gamma }^{22}}$ denote the reduced covariance matrices of subsystems A₁ and A₂, respectively; while ${{\Gamma }^{12}}$ and ${{\Gamma }^{21}}$ contain the expectation values of mixed quadratic terms. Let P₂ be the parity operator on subsystem A₂, and define the operators ${{\rho }_{+}}=\frac{1}{2}({{\rho }_{A}}+{{P}_{2}}{{\rho }_{A}}{{P}_{2}})$ and ${{\rho }_{-}}=\frac{1}{2}({{\rho }_{A}}-{{P}_{2}}{{\rho }_{A}}{{P}_{2}})$. By definition, we have that ${{\rho }_{A}}={{\rho }_{+}}+{{\rho }_{-}}$. Using the notation of section 3.1, ${{\rho }_{+}}$ and ${{\rho }_{-}}$ can be expanded as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\rho }_{+}} & = & \mathop{\sum }\limits_{\mathop{\underline{\kappa },\underline{\tau }}\limits_{|\underline{\tau }|{\rm even}}}{{w}_{\underline{\kappa },\underline{\tau }}}\;a_{{{m}_{1}}}^{{{\kappa }_{1}}}\;\ldots \;a_{{{m}_{2k}}}^{{{\kappa }_{2k}}}a_{{{n}_{1}}}^{{{\tau }_{1}}}\;\ldots \;a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}}, \\ {{\rho }_{-}} & = & \mathop{\sum }\limits_{\mathop{\underline{\kappa },\underline{\tau }}\limits_{|\underline{\tau }|\ {\rm odd}}}{{w}_{\underline{\kappa },\underline{\tau }}}\;a_{{{m}_{1}}}^{{{\kappa }_{1}}}\;\ldots \;a_{{{m}_{2k}}}^{{{\kappa }_{2k}}}a_{{{n}_{1}}}^{{{\tau }_{1}}}\;\ldots \;a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}}, \\ \end{array}\end{eqnarray} \tag{ 22 }$$

where the coefficients ${{w}_{\underline{\kappa },\underline{\tau }}}$ can be obtained from ${{\Gamma }_{A}}$ using the Wick rule, equation (5). By linearity of the partial transpose, $\rho _{A}^{{{T}_{2}}}=\rho _{+}^{{{T}_{2}}}+\rho _{-}^{{{T}_{2}}}$ follows, and $\rho _{\pm }^{{{T}_{2}}}$ can be obtained using equation (20):

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \rho _{+}^{{{T}_{2}}} & = & \mathop{\sum }\limits_{\mathop{\underline{\kappa },\underline{\tau }}\limits_{|\underline{\tau }|\ {\rm even}}}{{(-1)}^{|\underline{\tau }|/2}}{{w}_{\underline{\kappa },\underline{\tau }}}\;a_{{{m}_{1}}}^{{{\kappa }_{1}}}\ldots a_{{{m}_{2k}}}^{{{\kappa }_{2k}}}a_{{{n}_{1}}}^{{{\tau }_{1}}}\ldots a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}}, \\ \rho _{-}^{{{T}_{2}}} & = & \mathop{\sum }\limits_{\mathop{\underline{\kappa },\underline{\tau }}\limits_{|\underline{\tau }|\ {\rm odd}}}{{(-1)}^{(|\underline{\tau }|-1)/2}}{{w}_{\underline{\kappa },\underline{\tau }}}\;a_{{{m}_{1}}}^{{{\kappa }_{1}}}\ldots a_{{{m}_{2k}}}^{{{\kappa }_{2k}}}a_{{{n}_{1}}}^{{{\tau }_{1}}}\ldots a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}}. \\ \end{array}\end{eqnarray} \tag{ 23 }$$

Let us introduce the generalized Gaussian operators ${{O}_{+}}$ and ${{O}_{-}}$, with covariance matrices

$$\begin{eqnarray}{{\Gamma }_{+}}=\left( \begin{array}{ccccccccccccccc} {{\Gamma }^{11}} & {\rm i}{{\Gamma }^{12}} \\ {\rm i}{{\Gamma }^{21}} & -{{\Gamma }^{22}} \\ \end{array} \right),\qquad {{\Gamma }_{-}}=\left( \begin{array}{ccccccccccccccc} {{\Gamma }^{11}} & -{\rm i}{{\Gamma }^{12}} \\ -{\rm i}{{\Gamma }^{21}} & -{{\Gamma }^{22}} \\ \end{array} \right),\end{eqnarray} \tag{ 24 }$$

and consider the Majorana monomial expansion of these operators

$$\begin{eqnarray}&&{{O}_{\pm }}=\mathop{\sum }\limits_{\underline{\kappa },\underline{\tau }}o_{\underline{\kappa },\underline{\tau }}^{\pm }\;a_{{{m}_{1}}}^{{{\kappa }_{1}}}\ldots a_{{{m}_{2k}}}^{{{\kappa }_{2k}}}a_{{{n}_{1}}}^{{{\tau }_{1}}}\ldots a_{{{n}_{2\ell }}}^{{{\tau }_{2\ell }}}.\end{eqnarray} \tag{ 25 }$$

Since ${{O}_{+}}$ and ${{O}_{-}}$ are Gaussian operators, one can again obtain $o_{\underline{\kappa },\underline{\tau }}^{\pm }$ from ${{\Gamma }_{\pm }}$ using equation (5). Connecting the Wick-expansion form of ${{w}_{\underline{\kappa },\underline{\tau }}}$ with that of $o_{\underline{\kappa },\underline{\tau }}^{\pm }$, using the relation between ${{\Gamma }_{A}}$ and ${{\Gamma }_{\pm }}$, one can deduce that

$$\begin{eqnarray}o_{\underline{\kappa },\underline{\tau }}^{\pm }=\left\{ \begin{array}{ccccccccccccccc} \pm {\rm i}{{(-1)}^{(|\underline{\tau }|-1)/2}}{{w}_{\underline{\kappa },\underline{\tau }}} & {\rm when}|\underline{\tau }|{\rm odd}, \\ {{(-1)}^{|\underline{\tau }|/2}}{{w}_{\underline{\kappa },\underline{\tau }}} & {\rm when}|\underline{\tau }|{\rm even}. \\ \end{array} \right.\end{eqnarray} \tag{ 26 }$$

Comparing this with equation (23) it immediately follows that $\rho _{+}^{{{T}_{2}}}=\frac{1}{2}({{O}_{+}}+{{O}_{-}})$ and $\rho _{-}^{{{T}_{2}}}=\frac{{\rm i}}{2}({{O}_{-}}-{{O}_{+}})$. Thus, we obtain the decomposition

$$\begin{eqnarray}&&\rho _{A}^{{{T}_{2}}}=\frac{1-{\rm i}}{2}{{O}_{+}}+\frac{1+{\rm i}}{2}{{O}_{-}},\end{eqnarray} \tag{ 27 }$$

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

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, equation (27) can not be rewritten for the eigenvalues, and thus one does not have a direct access to the full spectrum of $\rho _{A}^{{{T}_{2}}}$. Nevertheless, there are a number of important properties which can be deduced by a direct investigation of the covariance matrices ${{\Gamma }_{\pm }}$. 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 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=\varnothing$ and a pure state on $A={{A}_{1}}\cup {{A}_{2}}$ given by $\rho =|\phi \rangle \langle \phi |$. Since for any pure Gaussian fermionic state ${{\Gamma }^{2}}=$ is satisfied, it follows that $[{{\Gamma }_{+}},{{\Gamma }_{-}}]$ = 0, which implies that ${{O}_{+}}$ and ${{O}_{-}}$ commute as well. Furthermore, since their eigenvalues are connected by a complex conjugation, the spectrum of ${{\rho }^{{{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 ${{\Gamma }_{+}}$. To obtain them, we first assume without loss of generality that $|{{A}_{1}}|\leqslant |{{A}_{2}}|$. Furthermore, for notational simplicity we also assume that $|A|$ is even, however, the results for the odd case follow trivially. Let $\pm {{\mu }_{k}}$ denote the eigenvalues of the reduced covariance matrix ${{\Gamma }^{11}}$ with $k=1,\ldots ,|{{A}_{1}}|$ and ${{\mu }_{k}}\geqslant 0$. As shown in appendix B, the eigenvalues $\pm \nu _{k}^{\pm }$ of ${{\Gamma }_{+}}$ can then be given as

$$\begin{eqnarray}\nu _{k}^{\pm }=\left\{ \begin{array}{ccccccccccccccc} {{\mu }_{k}}\pm {\rm i}\sqrt{1-\mu _{k}^{2}} & {\rm when}\ k=1,\ldots ,|{{A}_{1}}|, \\ 1 & {\rm when}\ k=|{{A}_{1}}|+1,\ldots ,|A|/2. \\ \end{array} \right.\end{eqnarray} \tag{ 28 }$$

The canonical diagonalized form of the Gaussian operator ${{O}_{+}}$ reads

$$\begin{eqnarray}&&{{O}_{+}}=\mathop{\prod }\limits_{k=1}^{|A|/2}\mathop{\prod }\limits_{{{\sigma }_{k}}=\pm }\frac{+{\rm i}\nu _{k}^{{{\sigma }_{k}}}b_{2k-1}^{{{\sigma }_{k}}}b_{k}^{{{\sigma }_{k}}}}{2},\end{eqnarray} \tag{ 29 }$$

where $b_{j}^{\pm }$ are Majorana operators obtained from a_j via the orthogonal transformation which diagonalizes ${{\Gamma }_{+}}$. 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,\ldots ,|{{A}_{1}}|$ as

$$\begin{eqnarray}&&\omega _{k}^{{{\sigma }_{k}}\sigma _{k}^{\prime }}=\frac{1+{{\sigma }_{k}}\nu _{k}^{+}}{2}\frac{1+\sigma _{k}^{\prime }\nu _{k}^{-}}{2},\end{eqnarray} \tag{ 30 }$$

with ${{\sigma }_{k}}=\pm$ and $\sigma _{k}^{\prime }=\pm$. Using equation (28) this yields

$$\begin{eqnarray}&&\omega _{k}^{++}=\frac{1+{{\mu }_{k}}}{2},\quad \omega _{k}^{--}=\frac{1-{{\mu }_{k}}}{2},\quad \omega _{k}^{+-}=-\omega _{k}^{-+}=\frac{{\rm i}}{2}\sqrt{1-\mu _{k}^{2}},\end{eqnarray} \tag{ 31 }$$

where we used the property $\nu _{k}^{+}\nu _{k}^{-}=1$. The nonzero eigenvalues of ${{O}_{+}}$ can then be written down as

$$\begin{eqnarray}&&{{\Omega }_{\underline{\sigma }\,{{{\underline{\sigma }}}^{\prime }}}}=\mathop{\prod }\limits_{k=1}^{|{{A}_{1}}|}\omega _{k}^{{{\sigma }_{k}}\sigma _{k}^{\prime }}=\mathop{\prod }\limits_{{{\sigma }_{k}}=\sigma _{k}^{\prime }}\frac{1+{{\sigma }_{k}}{{\mu }_{k}}}{2}\mathop{\prod }\limits_{{{\sigma }_{k}}=-\sigma _{k}^{\prime }}\frac{{{\sigma }_{k}}{\rm i}\sqrt{1-\mu _{k}^{2}}}{2},\end{eqnarray} \tag{ 32 }$$

where $\underline{\sigma }$ and $\underline{\sigma }^{\prime}$ are the signature arrays corresponding to the eigenvalue. Note that the additional $\nu _{k}^{\pm }=1$ in equation (28) lead to further eigenvalues of ${{O}_{+}}$ that are all equal to zero.

The products in equation (32) are either real or purely imaginary and the eigenvalues of ${{\rho }^{{{T}_{2}}}}$ thus follow as $\operatorname{Re}{{\Omega }_{\underline{\sigma }{{{\underline{\sigma }}}^{\prime }}}}$ or ${\rm Im}{{\Omega }_{\underline{\sigma }{{{\underline{\sigma }}}^{\prime }}}}$, respectively. It is instructive, however, to derive the same spectrum using the Schmidt decomposition of $|\phi \rangle$. Dividing the Hilbert space $\mathcal{H}={{\mathcal{H}}_{1}}\otimes {{\mathcal{H}}_{2}}$ into two parts⁴ , one has

$$\begin{eqnarray}&&|\phi \rangle =\mathop{\sum }\limits_{i}\sqrt{{{\lambda }_{i}}}|\phi _{i}^{1}\rangle |\phi _{i}^{2}\rangle ,\end{eqnarray} \tag{ 33 }$$

with the RDM eigenvalues ${{\lambda }_{i}}$ of subsystem A₁. Clearly, the state is supported on a smaller Hilbert space ${{\mathcal{H}}_{1}}\otimes {{\mathcal{H}}_{1}}$ and is invariant under the action of a flip operation defined by $|\phi _{i}^{1}\rangle |\phi _{j}^{2}\rangle \to |\phi _{j}^{1}\rangle |\phi _{i}^{2}\rangle$. The partial transpose of ρ,

$$\begin{eqnarray}&&{{\rho }^{{{T}_{2}}}}=(|\phi \rangle {{\langle \phi |)}^{{{T}_{2}}}}=\mathop{\sum }\limits_{i,j}\sqrt{{{\lambda }_{i}}{{\lambda }_{j}}}|\phi _{i}^{1}\rangle \langle \phi _{j}^{1}|\otimes |\phi _{j}^{2}\rangle \langle \phi _{i}^{2}|,\end{eqnarray} \tag{ 34 }$$

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

$$\begin{eqnarray}&&{{\rho }^{{{T}_{2}}}}|\phi _{i}^{1}\rangle |\phi _{i}^{2}\rangle ={{\lambda }_{i}}|\phi _{i}^{1}\rangle |\phi _{i}^{2}\rangle ,\qquad {{\rho }^{{{T}_{2}}}}|\phi _{ij}^{\pm }\rangle =\pm \sqrt{{{\lambda }_{i}}{{\lambda }_{j}}}|\phi _{ij}^{\pm }\rangle ,\end{eqnarray} \tag{ 35 }$$

where we introduced the notation

$$\begin{eqnarray}&&|\phi _{ij}^{\pm }\rangle =\frac{1}{\sqrt{2}}(|\phi _{i}^{1}\rangle |\phi _{j}^{2}\rangle \pm |\phi _{j}^{1}\rangle |\phi _{i}^{2}\rangle )\qquad i\ne j.\end{eqnarray} \tag{ 36 }$$

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

$$\begin{eqnarray}&&\sqrt{{{\lambda }_{{\underline{\sigma }}}}{{\lambda }_{{{{\underline{\sigma }}}^{\prime }}}}}=\mathop{\prod }\limits_{{{\sigma }_{k}}=\sigma _{k}^{\prime }}\frac{1+{{\sigma }_{k}}{{\mu }_{k}}}{2}\mathop{\prod }\limits_{{{\sigma }_{k}}=-\sigma _{k}^{\prime }}\frac{\sqrt{1-\mu _{k}^{2}}}{2},\end{eqnarray} \tag{ 37 }$$

where the signature $\underline{\sigma }$ has again components ${{\sigma }_{k}}=\pm$. Comparing equation (37) to (32), one indeed recognizes ${{\Omega }_{\underline{\sigma }{{{\underline{\sigma }}}^{\prime }}}}$ 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

$$\begin{eqnarray}&&{\rm Tr}|{{\rho }^{{{T}_{2}}}}|=1-2{\rm T}{{{\rm r}}_{{\rm o}}}\;{{\rho }^{{{T}_{2}}}}=\mathop{\prod }\limits_{k}\left( 1+\sqrt{1-\mu _{k}^{2}} \right),\end{eqnarray} \tag{ 38 }$$

where we introduced the notation ${\rm T}{{{\rm r}}_{{\rm e}/{\rm o}}}$ for traces taken over the even/odd subspace. Thus, the pure-state logarithmic negativity is given by the simple formula

$$\begin{eqnarray}&&\mathcal{E}={\rm ln} \;{\rm Tr}|{{\rho }^{{{T}_{2}}}}|=\mathop{\sum }\limits_{k}{\rm ln} \left( 1+\sqrt{1-\mu _{k}^{2}} \right).\end{eqnarray} \tag{ 39 }$$

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

$$\begin{eqnarray}&&{{S}_{n}}=\frac{1}{1-n}\mathop{\sum }\limits_{k}{\rm ln} \left[ {{\left( \frac{1+{{\mu }_{k}}}{2} \right)}^{n}}+\ {{\left( \frac{1-{{\mu }_{k}}}{2} \right)}^{n}} \right],\end{eqnarray} \tag{ 40 }$$

the well-known equality $\mathcal{E}={{S}_{1/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 $|{{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, 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 ${{\Gamma }_{+}}$ is invariant with respect to a sign change and complex conjugation, hence the eigenvalues can be collected into two families of quadruplets

$$\begin{eqnarray}&&\left\{ {{z}_{k}},z_{k}^{*},-{{z}_{k}},-z_{k}^{*} \right\},\quad k\in ({\rm I})\qquad \left\{ {\rm i}{{u}_{k}},-{\rm i}{{v}_{k}},-{\rm i}{{u}_{k}},{\rm i}{{v}_{k}} \right\},\quad k\in ({\rm II}),\end{eqnarray} \tag{ 41 }$$

where in family (I) we choose $\operatorname{Re}{{z}_{k}}\gt 0$ and ${\rm Im}{{z}_{k}}\gt 0$, whereas ${{u}_{k}}\gt {{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}}\ne {{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\lt 1$ in equation (7).

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

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \omega _{k}^{{{\sigma }_{k}}\sigma _{k}^{\prime }} & = & \left\{ \begin{array}{ccccccccccccccc} \frac{1+|{{z}_{k}}{{|}^{2}}+{{\sigma }_{k}}2\operatorname{Re}{{z}_{k}}}{4}, & {{\sigma }_{k}}=\sigma _{k}^{\prime } \\ \frac{1-|{{z}_{k}}{{|}^{2}}+{{\sigma }_{k}}2{\rm iIm}{{z}_{k}}}{4}, & {{\sigma }_{k}}=-\sigma _{k}^{\prime } \\ \end{array} \right.,\quad k\in ({\rm I}) \\ \omega _{k}^{{{\sigma }_{k}}\sigma _{k}^{\prime }} & = & \left\{ \begin{array}{ccccccccccccccc} \frac{1+{{u}_{k}}{{v}_{k}}+{{\sigma }_{k}}{\rm i}(u-v)}{4}, & {{\sigma }_{k}}=\sigma _{k}^{\prime } \\ \frac{1-{{u}_{k}}{{v}_{k}}+{{\sigma }_{k}}{\rm i}(u+v)}{4}, & {{\sigma }_{k}}=-\sigma _{k}^{\prime } \\ \end{array} \right.,\quad k\in ({\rm II}) \\ \end{array}\end{eqnarray} \tag{ 42 }$$

within each quadruplet, and the eigenvalues ${{\Omega }_{\underline{\sigma }\,{{{\underline{\sigma }}}^{\prime }}}}$ of ${{O}_{+}}$ are again given in the factorized form of equation (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 ${{\rho }^{{{T}_{2}}}}$. 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 $R\sigma _{n}^{\alpha }{{R}^{\dagger }}=\sigma _{N-n}^{\alpha }$ (with $\alpha =x,y,z$), implying the action $Rc_{n}^{\dagger }{{R}^{\dagger }}=Pc_{N-n}^{\dagger }$ on the creation operators, where P is the parity operator. Using this, it follows that the sign factor associated to an eigenvector of ${{O}_{+}}$ reads

$$\begin{eqnarray}{{S}_{\underline{\sigma }\,{{{\underline{\sigma }}}^{\prime }}}}=\operatorname{Re}\mathop{\prod }\limits_{k}{{s}_{{{\sigma }_{k}}\sigma _{k}^{\prime }}}+{\rm Im}\mathop{\prod }\limits_{k}{{s}_{{{\sigma }_{k}}\sigma _{k}^{\prime }}},\qquad {{s}_{{{\sigma }_{k}}\sigma _{k}^{\prime }}}=\left\{ \begin{array}{ccccccccccccccc} 1, & {\rm if}\ {{\sigma }_{k}}=\sigma _{k}^{\prime }, \\ {{\sigma }_{k}}i, & {\rm if}\ {{\sigma }_{k}}=-\sigma _{k}^{\prime }, \\ \end{array} \right.\end{eqnarray} \tag{ 43 }$$

which can also be verified by considering the pure state limit. The parity of the ${{O}_{-}}$ eigenvectors are simply obtained through the factors $s_{{{\sigma }_{k}}\sigma _{k}^{\prime }}^{*}$.

With the knowledge of ${{S}_{\underline{\sigma }\,{{{\underline{\sigma }}}^{\prime }}}}$, we are now able to carry out signed traces of the form

$$\begin{eqnarray}&&{\rm T}{{{\rm r}}_{{\rm e}}}{{O}_{+}}-{\rm T}{{{\rm r}}_{{\rm o}}}{{O}_{+}}=\mathop{\sum }\limits_{\underline{\sigma },{{{\underline{\sigma }}}^{\prime }}}{{S}_{\underline{\sigma }\,{{{\underline{\sigma }}}^{\prime }}}}{{\Omega }_{\underline{\sigma }\,{{{\underline{\sigma }}}^{\prime }}}}.\end{eqnarray} \tag{ 44 }$$

Note that the terms in the sum of equation (44) completely factorize in the quadruplet index k. Furthermore, using the fact that ${\rm T}{{{\rm r}}_{{\rm e}/{\rm o}}}{{O}_{-}}={{({\rm T}{{{\rm r}}_{{\rm e}/{\rm o}}}{{O}_{+}})}^{*}}$ and thus ${\rm T}{{{\rm r}}_{{\rm e}/{\rm o}}}{{\rho }^{{{T}_{2}}}}={\rm ReT}{{{\rm r}}_{{\rm e}/{\rm o}}}{{O}_{+}}+{\rm ImT}{{{\rm r}}_{{\rm e}/{\rm o}}}{{O}_{+}}$, a simple calculation leads to

$$\begin{eqnarray}&&1-2{\rm T}{{{\rm r}}_{{\rm o}}}{{\rho }^{{{T}_{2}}}}={\rm T}{{{\rm r}}_{{\rm e}}}{{\rho }^{{{T}_{2}}}}-{\rm T}{{{\rm r}}_{{\rm o}}}{{\rho }^{{{T}_{2}}}}=\mathop{\prod }\limits_{k\in ({\rm I})}\frac{1+|{{z}_{k}}{{|}^{2}}+2{\rm Im}{{z}_{k}}}{2}\mathop{\prod }\limits_{k\in ({\rm II})}\frac{1+{{u}_{k}}{{v}_{k}}+{{u}_{k}}+{{v}_{k}}}{2}.\end{eqnarray} \tag{ 45 }$$

Finally, we define the quantity

$$\begin{eqnarray}&&{{\mathcal{E}}_{{\rm o}}}={\rm ln} \;{\rm max} (1-2{\rm T}{{{\rm r}}_{{\rm o}}}{{\rho }^{{{T}_{2}}}},1),\end{eqnarray} \tag{ 46 }$$

which clearly gives a lower bound for the logarithmic negativity, $\mathcal{E}\geqslant {{\mathcal{E}}_{{\rm 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 [50, 51].

To test the bound ${{\mathcal{E}}_{{\rm o}}}$ against the exact value of the logarithmic negativity, we considered small TI chains with $N\leqslant 10$ and obtained $\mathcal{E}$ via exact diagonalization of ${{\rho }^{{{T}_{2}}}}$. 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, ${{\mathcal{E}}_{{\rm o}}}$ indeed exactly coincides with $\mathcal{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 $\mathcal{E}-{{\mathcal{E}}_{{\rm o}}}$. Nevertheless, for the small system sizes considered, it appears that ${{\mathcal{E}}_{{\rm o}}}$ gives a very good approximation up to temperatures $T\approx 1/N$, above which it starts to deviate significantly.

Zoom In Zoom Out Reset image size

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 ${{\rho }_{A}}$ after tracing out the sites of B is a mixed Gaussian state, associated to the reduced covariance matrix ${{\Gamma }_{A}}$, with indices running over sites in A [52].

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]

$$\begin{eqnarray}&&\mathcal{E}(\ell ,L)=\frac{c}{4}{\rm ln} \left[ \frac{L}{\pi }\frac{{{{\rm sin} }^{2}}\left( \frac{\pi \ell }{L} \right)}{{\rm sin} \left( \frac{2\pi \ell }{L} \right)} \right]+{\rm const}.,\end{eqnarray} \tag{ 47 }$$

with the central charge c and a non-universal constant. However, subtracting the value at $\ell =L/4$ one obtains a universal scaling function

$$\begin{eqnarray}&&\epsilon (z)=\mathcal{E}(\ell ,L)-\mathcal{E}(L/4,L)=\frac{1}{8}{\rm ln} [{\rm tan} (\pi z)],\end{eqnarray} \tag{ 48 }$$

where $z=\ell /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 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 ${{\mathcal{E}}_{{\rm o}}}$ is given by the same formula in terms of the eigenvalues of ${{\Gamma }_{+}}$. 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 $\rho _{A}^{{{T}_{2}}}$ 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 $\rho _{A}^{{{T}_{2}}}$ 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 ${{\mathcal{E}}_{{\rm o}}}$ is tight.

In spite of the above findings, ${{\mathcal{E}}_{{\rm o}}}$ 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 ${{\mathcal{E}}_{{\rm o}}}(L/4,L)$ against h on the left of figure 3. Furthermore, defining the quantity ${{\epsilon }_{{\rm o}}}$ analogously to equation (48), one finds an excellent data collapse when plotted against the variable $z=\ell /L$, see right of figure 3. The scaling function is found to be given by

$$\begin{eqnarray}&&{{\epsilon }_{{\rm o}}}(z)={{\mathcal{E}}_{{\rm o}}}(\ell ,L)-{{\mathcal{E}}_{{\rm o}}}(L/4,L)=\frac{1}{16}{\rm ln} \left[ \frac{{\rm tan} (\pi z)}{2{{{\rm cos} }^{2}}(\pi z)} \right].\end{eqnarray} \tag{ 49 }$$

Zoom In Zoom Out Reset image size

Although the functional form of ${{\epsilon }_{{\rm 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 $\epsilon (z)$ in equation (48), however, the argument is modified as well. We also performed a calculation directly in the thermodynamic limit $L\to \infty$, and found ${{\mathcal{E}}_{{\rm o}}}(\ell )=1/16{\rm ln} \ell +{\rm 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}}=\frac{1}{2}({{\delta }_{m,n-1}}+{{\delta }_{m,n+1}})$. In complete analogy with the result for the bipartite entanglement [53], one finds $\mathcal{E}_{{\rm o}}^{XX}(2\ell )=2\mathcal{E}_{{\rm o}}^{{\rm TI}}(\ell )$ for h = 1, and thus a doubled prefactor $1/8$ with respect to the critical TI chain. Therefore, even though ${{\mathcal{E}}_{{\rm o}}}$ does not approximate $\mathcal{E}$ well, it shows exactly the same universal behaviour, suggesting it as an entanglement indicator which is extremely simple to calculate.

Our result for the partial transpose, equation (27), can be further tested by looking at traces of integer powers of the partial transpose $\rho _{A}^{{{T}_{2}}}$ which can also be carried out within CFT. Since the identity ${\rm Tr}{{(\rho _{A}^{{{T}_{2}}})}^{2}}={\rm Tr}\rho _{A}^{2}$ 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]

$$\begin{eqnarray}&&{{R}_{3}}(\ell ,L)={\rm lnTr}{{(\rho _{A}^{{{T}_{2}}})}^{3}}=-\frac{1}{9}{\rm ln} \left[ \frac{{{L}^{3}}}{{{\pi }^{3}}}{{{\rm sin} }^{2}}\left( \frac{\pi \ell }{L} \right){\rm sin} \left( \frac{2\pi \ell }{L} \right) \right]+{\rm const}.\end{eqnarray} \tag{ 50 }$$

Similarly to equation (48), a universal scaling function can be defined as [24]

$$\begin{eqnarray}&&{{r}_{3}}(z)={{R}_{3}}(\ell ,L)-{{R}_{3}}(L/4,L)=-\frac{1}{9}{\rm ln} \left[ 2{{{\rm sin} }^{2}}(\pi z){\rm sin} (2\pi z) \right],\end{eqnarray} \tag{ 51 }$$

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

$$\begin{eqnarray}&&{{R}_{3}}(\ell ,\beta )=-\frac{1}{9}{\rm ln} \left[ \frac{{{\beta }^{3}}}{{{\pi }^{3}}}{{{\rm sinh} }^{2}}\left( \frac{\pi \ell }{\beta } \right){\rm sinh} \left( \frac{2\pi \ell }{\beta } \right) \right]+{\rm const}.\end{eqnarray} \tag{ 52 }$$

We now show how the above traces can be calculated with the covariance matrix formalism. Expanding the third power of $\rho _{A}^{{{T}_{2}}}$ in equation (27) and taking the trace one arrives at

$$\begin{eqnarray}&&{\rm Tr}{{(\rho _{A}^{{{T}_{2}}})}^{3}}=-\frac{1}{2}{\rm Tr}(O_{+}^{3})+\frac{3}{2}{\rm Tr}(O_{+}^{2}{{O}_{-}}),\end{eqnarray} \tag{ 53 }$$

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

$$\begin{eqnarray}&&{\rm Tr}{{(\rho _{A}^{{{T}_{2}}})}^{3}}=\mp \frac{1}{2}\sqrt{{\rm det} \left( \frac{1+3\Gamma _{+}^{2}}{4} \right)}+\frac{3}{2}\sqrt{{\rm det} \left( \frac{1+\Gamma _{+}^{2}+2{{\Gamma }_{+}}{{\Gamma }_{-}}}{4} \right)},\end{eqnarray} \tag{ 54 }$$

where the sign of the first term depends on the spectrum of ${{\Gamma }_{+}}$, 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 ${{\Gamma }_{\pm }}$ 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.

Zoom In Zoom Out Reset image size

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 ${{\Gamma }_{\pm }}$. 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 $\rho _{A}^{{{T}_{2}}}$ 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 ${{\mathcal{E}}_{{\rm 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 of the traces for such a case is still realizable, but presumably more tedious.

It would also be interesting to see whether the lower bound ${{\mathcal{E}}_{{\rm 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 ${{\mathcal{E}}_{{\rm o}}}$ and $\mathcal{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.

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 $[{{\Gamma }_{+}},{{\Gamma }_{-}}]=0$ is satisfied, an example being the isotropic Gaussian states [44], for which the covariance matrix satisfies ${{\Gamma }^{2}}={{\lambda }^{2}}$ with some $0\leqslant \lambda \lt 1$. The situation is then analogous to the pure-state case and the corresponding calculation can be generalized, which we leave for future studies.

We would like to thank Leonardo Banchi and Fernando Brandão 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).

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

$$\begin{eqnarray}&&\rho =\mathop{\prod }\limits_{k=1,2}\frac{+{\rm i}{{\nu }_{k}}{{b}_{2k-1}}{{b}_{2k}}}{2},\end{eqnarray} \tag{ A.1 }$$

where ${{\nu }_{k}}\in [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 equation (13), and assume also reflection symmetry. These states are parametrized by a single angle θ beside the covariance matrix eigenvalues ${{\nu }_{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

$$\begin{eqnarray}{{\rho }^{{{T}_{2}}}}=\frac{}{4}+\frac{{{\nu }_{1}}{{\nu }_{2}}}{4}\sigma _{1}^{z}\sigma _{2}^{z}+\frac{{{\nu }_{1}}+{{\nu }_{2}}}{4}\left( \begin{array}{ccccccccccccccc} {\rm cos} 2\theta & {} & {} & 0 \\ {} & 0 & {\rm sin} 2\theta & {} \\ {} & {\rm sin} 2\theta & 0 & {} \\ 0 & {} & {} & -{\rm cos} 2\theta \\ \end{array} \right)+\frac{{{\nu }_{1}}-{{\nu }_{2}}}{4}\left( \begin{array}{ccccccccccccccc} 0 & {} & {} & 1 \\ {} & 0 & 0 & {} \\ {} & 0 & 0 & {} \\ 1 & {} & {} & 0 \\ \end{array} \right).\end{eqnarray} \tag{ A.2 }$$

Note that, besides the diagonals, all matrix elements vanish and are thus not shown. It is straightforward to obtain the four eigenvalues

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\lambda }_{1,2}} & = & \left[ 1+{{\nu }_{1}}{{\nu }_{2}}\pm \sqrt{{{({{\nu }_{1}}+{{\nu }_{2}})}^{2}}{{{\rm cos} }^{2}}2\theta +{{({{\nu }_{1}}-{{\nu }_{2}})}^{2}}} \right]/4, \\ {{\lambda }_{3,4}} & = & \left[ 1-{{\nu }_{1}}{{\nu }_{2}}\pm ({{\nu }_{1}}+{{\nu }_{2}}){\rm sin} 2\theta \right]/4. \\ \end{array}\end{eqnarray} \tag{ A.3 }$$

Using the parity operator ${{P}_{2}}=\sigma _{2}^{z}$, one can also construct $\rho _{\pm }^{{{T}_{2}}}=({{\rho }^{{{T}_{2}}}}\pm {{P}_{2}}{{\rho }^{{{T}_{2}}}}{{P}_{2}})/2$ with matrix elements

$$\begin{eqnarray}\rho _{+}^{{{T}_{2}}}=\frac{}{4}+\frac{{{\nu }_{1}}{{\nu }_{2}}}{4}\sigma _{1}^{z}\sigma _{2}^{z}+\frac{{{\nu }_{1}}+{{\nu }_{2}}}{4}\left( \begin{array}{ccccccccccccccc} {\rm cos} 2\theta & {} & {} & 0 \\ {} & 0 & 0 & {} \\ {} & 0 & 0 & {} \\ 0 & {} & {} & -{\rm cos} 2\theta \\ \end{array} \right),\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}\rho _{-}^{{{T}_{2}}}=\frac{{{\nu }_{1}}+{{\nu }_{2}}}{4}\left( \begin{array}{ccccccccccccccc} 0 & {} & {} & 0 \\ {} & 0 & {\rm sin} 2\theta & {} \\ {} & {\rm sin} 2\theta & 0 & {} \\ 0 & {} & {} & 0 \\ \end{array} \right)+\frac{{{\nu }_{1}}-{{\nu }_{2}}}{4}\left( \begin{array}{ccccccccccccccc} 0 & {} & {} & 1 \\ {} & 0 & 0 & {} \\ {} & 0 & 0 & {} \\ 1 & {} & {} & 0 \\ \end{array} \right).\end{eqnarray} \tag{ A.5 }$$

The eigenvalues of the operator $\rho _{+}^{{{T}_{2}}}+{\rm i}\rho _{-}^{{{T}_{2}}}$ then read

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\Omega }_{1,2}} & = & \left[ 1+{{\nu }_{1}}{{\nu }_{2}}\pm \sqrt{{{({{\nu }_{1}}+{{\nu }_{2}})}^{2}}{{{\rm cos} }^{2}}2\theta -{{({{\nu }_{1}}-{{\nu }_{2}})}^{2}}} \right]/4, \\ {{\Omega }_{3,4}} & = & \left[ 1-{{\nu }_{1}}{{\nu }_{2}}\pm {\rm i}({{\nu }_{1}}+{{\nu }_{2}}){\rm sin} 2\theta \right]/4. \\ \end{array}\end{eqnarray} \tag{ A.6 }$$

Note that we have ${{\lambda }_{3,4}}=\operatorname{Re}{{\Omega }_{3,4}}+{\rm Im}{{\Omega }_{3,4}}$ which simply follows from the fact that $\rho _{+}^{{{T}_{2}}}$ and $\rho _{-}^{{{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 ${{\Gamma }_{+}}$ has indeed eigenvalues given by equation (A.6). The covariance matrix Γ for the Gaussian state (A.1) and the associated ${{\Gamma }_{+}}$ have the form

$$\begin{eqnarray}\Gamma ={\rm i}\left( \begin{array}{ccccccccccccccc} 0 & c & 0 & {{s}_{-}} \\ -c & 0 & {{s}_{+}} & 0 \\ 0 & -{{s}_{+}} & 0 & c \\ -{{s}_{-}} & 0 & -c & 0 \\ \end{array} \right),\qquad {{\Gamma }_{+}}={\rm i}\left( \begin{array}{ccccccccccccccc} 0 & c & 0 & {\rm i}{{s}_{-}} \\ -c & 0 & {\rm i}{{s}_{+}} & 0 \\ 0 & -{\rm i}{{s}_{+}} & 0 & -c \\ -{\rm i}{{s}_{-}} & 0 & c & 0 \\ \end{array} \right),\end{eqnarray} \tag{ A.7 }$$

with the shorthand notation

$$\begin{eqnarray}&&c=\frac{{{\nu }_{1}}+{{\nu }_{2}}}{2}{\rm cos} 2\theta ,\qquad {{s}_{\pm }}=\frac{{{\nu }_{1}}+{{\nu }_{2}}}{2}{\rm sin} 2\theta \pm \frac{{{\nu }_{1}}-{{\nu }_{2}}}{2}.\end{eqnarray} \tag{ A.8 }$$

The four eigenvalues $\pm {{\nu }^{\pm }}$ of ${{\Gamma }_{+}}$ can be computed with

$$\begin{eqnarray}&&{{\nu }^{\pm }}=\sqrt{{{c}^{2}}-{{\left( \frac{{{s}_{+}}-{{s}_{-}}}{2} \right)}^{2}}}\pm {\rm i}\frac{{{s}_{+}}+{{s}_{-}}}{2}.\end{eqnarray} \tag{ A.9 }$$

If the operator ${{O}_{+}}$ is Gaussian, its eigenvalues must have the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\Omega }_{++}} & = & \frac{1+{{\nu }_{+}}}{2}\frac{1+{{\nu }_{-}}}{2},\quad {{\Omega }_{--}}=\frac{1-{{\nu }_{+}}}{2}\frac{1-{{\nu }_{-}}}{2}, \\ {{\Omega }_{+-}} & = & \frac{1+{{\nu }_{+}}}{2}\frac{1-{{\nu }_{-}}}{2},\quad {{\Omega }_{-+}}=\frac{1-{{\nu }_{+}}}{2}\frac{1+{{\nu }_{-}}}{2}. \\ \end{array}\end{eqnarray} \tag{ A.10 }$$

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 ${{\rho }^{{{T}_{2}}}}$. comparing the expectation values ${\rm Tr}({{\rho }^{{{T}_{2}}}}{{a}_{m}}{{a}_{n}})$, where $m,n=1,\ldots ,4$, with ${\rm Tr}({{\rho }^{{{T}_{2}}}}{{a}_{1}}{{a}_{2}}{{a}_{3}}{{a}_{4}})$, one observes that the Wick expansion, equation (5), does not hold, unless ${{c}^{2}}={{\nu }_{1}}{{\nu }_{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.

Here we consider the eigenvalue problem of the modified covariance matrices ${{\Gamma }_{\pm }}$, that are associated to a pure-state covariance matrix Γ of the form (13). The RDMs for subsystems A₁ and A₂ are determined via the reduced covariance matrices ${{\Gamma }^{11}}$ and ${{\Gamma }^{22}}$. Since they split into two submatrices, one could equivalently consider the squared eigenvalue problem of the matrices

$$\begin{eqnarray}G_{mn}^{\alpha \alpha }=\mathop{\sum }\limits_{l\in \alpha }{{g}_{ml}}{{g}_{nl}}=\mathop{\sum }\limits_{p,q}M_{pq}^{\alpha }\psi _{p}^{\alpha }(m)\psi _{q}^{\alpha }(n),\qquad \psi _{q}^{\alpha }(m)=\left\{ \begin{array}{ccccccccccccccc} {{\psi }_{q}}(m) & m\in {{A}_{\alpha }}, \\ 0 & m\notin {{A}_{\alpha }}, \\ \end{array} \right.\end{eqnarray} \tag{ B.1 }$$

with nonzero matrix elements only within $m,n\in {{A}_{\alpha }}$, $\alpha =1,2$. The overlap matrices ${{M}^{\alpha }}$ have matrix elements

$$\begin{eqnarray}M_{pq}^{\alpha }=\mathop{\sum }\limits_{l}\phi _{p}^{\alpha }(l)\phi _{q}^{\alpha }(l),\qquad \phi _{q}^{\alpha }(l)=\left\{ \begin{array}{ccccccccccccccc} {{\phi }_{q}}(l) & l\in {{A}_{\alpha }}, \\ 0 & l\notin {{A}_{\alpha }}. \\ \end{array} \right.\end{eqnarray} \tag{ B.2 }$$

Note that both G¹¹ and G²² have the same eigenvalues $\mu _{k}^{2}\leqslant 1$ with $k=1,...,{\rm min} (|{{A}_{1}}|,|{{A}_{2}}|)$, whereas $\mu _{k}^{2}=1$ for the remaining eigenvalues of the larger matrix. We also introduce the block-diagonal matrix

$$\begin{eqnarray}G=\left( \begin{array}{ccccccccccccccc} {{G}^{11}} & 0 \\ 0 & {{G}^{22}} \\ \end{array} \right),\end{eqnarray} \tag{ B.3 }$$

with all the nontrivial eigenvalues being doubly degenerate.

The matrix elements of the covariance matrices ${{\Gamma }_{\pm }}$ in equation (24) are determined through

$$\begin{eqnarray}&&g_{mn}^{\pm }=\mathop{\sum }\limits_{q}\psi _{q}^{\pm }(m)\phi _{q}^{\pm }(n),\end{eqnarray} \tag{ B.4 }$$

with the vectors

$$\begin{eqnarray}&&\phi _{q}^{\pm }(l)=\phi _{q}^{1}(l)\pm {\rm i}\phi _{q}^{2}(l),\qquad \psi _{q}^{\pm }(l)=\psi _{q}^{1}(l)\pm {\rm i}\psi _{q}^{2}(l).\end{eqnarray} \tag{ B.5 }$$

Thus the spectrum of ${{\Gamma }_{\pm }}$ follows from the eigenvalues of the squared matrix

$$\begin{eqnarray}&&G_{mn}^{\pm }=\mathop{\sum }\limits_{p,q}(M_{pq}^{1}-M_{pq}^{2})\psi _{p}^{\pm }(m)\psi _{q}^{\pm }(n).\end{eqnarray} \tag{ B.6 }$$

Inserting (B.5) and using the completeness property $Mpq+M_{pq}^{2}={{\delta }_{pq}}$, the matrices ${{G}^{\pm }}$ have the block form

$$\begin{eqnarray}{{G}^{\pm }}=\left( \begin{array}{ccccccccccccccc} 2{{G}^{11}}- & \pm 2{\rm i}F \\ \pm 2{\rm i}{{F}^{T}} & 2{{G}^{22}}- \\ \end{array} \right),\end{eqnarray} \tag{ B.7 }$$

with

$$\begin{eqnarray}&&{{F}_{mn}}=\mathop{\sum }\limits_{p,q}M_{pq}^{1}\psi _{p}^{1}(m)\psi _{q}^{2}(n).\end{eqnarray} \tag{ B.8 }$$

It is easy to check that F satisfies

$$\begin{eqnarray}&&F{{F}^{T}}={{G}^{11}}(-{{G}^{11}}),\qquad {{F}^{T}}F={{G}^{22}}(-{{G}^{22}}),\end{eqnarray} \tag{ B.9 }$$

and thus the following matrix identity holds

$$\begin{eqnarray}&&{{({{G}^{\pm }}-2\;G+)}^{2}}=-4G(-G).\end{eqnarray} \tag{ B.10 }$$

Rewriting in terms of the eigenvalues ${{(\nu _{k}^{\pm })}^{2}}$ and $\mu _{k}^{2}$ of ${{G}^{\pm }}$ and G, respectively, one arrives at

$$\begin{eqnarray}&&{{(\nu _{k}^{\pm })}^{2}}=2\mu _{k}^{2}-1\pm 2{\rm i}{{\mu }_{k}}\sqrt{1-\mu _{k}^{2}}={{\left( {{\mu }_{k}}\pm {\rm i}\sqrt{1-\mu _{k}^{2}} \right)}^{2}}.\end{eqnarray} \tag{ B.11 }$$

Let us consider the Gaussian operators ${{O}_{\pm }}$ corresponding to the generalized covariance matrices ${{\Gamma }_{\pm }}={\rm tanh} ({{W}_{\pm }}/2)$. With a denoting the vector of Majorana operators, we introduce

$$\begin{eqnarray}&&O\left( {{W}_{\pm }} \right)={\rm exp} \left( \frac{a{{W}_{\pm }}a}{4} \right),\qquad Z({{W}_{\pm }})={\rm Tr}\;{\rm exp} \left( \frac{a{{W}_{\pm }}a}{4} \right),\end{eqnarray} \tag{ C.1 }$$

and thus ${{O}_{\pm }}=O\left( {{W}_{\pm }} \right)/Z({{W}_{\pm }})$ Our aim is to calculate various traces of the form ${\rm Tr}(O_{+}^{m}O_{-}^{n})$ with some integers m and n. Following the lines of [48], we first introduce the notation

$$\begin{eqnarray}&&\left\{ {{W}_{1}},{{W}_{2}} \right\}={\rm Tr}\left[ O\left( {{W}_{1}} \right)O\left( {{W}_{2}} \right) \right].\end{eqnarray} \tag{ C.2 }$$

We also note the simple fact that ${{O}^{m}}\left( {{W}_{\pm }} \right)=O\left( m{{W}_{\pm }} \right)$. Hence, the traces we consider can be written in the form

$$\begin{eqnarray}&&{\rm Tr}(O_{+}^{m})=\frac{Z(m{{W}_{+}})}{{{Z}^{m}}({{W}_{+}})},\qquad {\rm Tr}(O_{+}^{m}O_{-}^{n})=\frac{\left\{ m{{W}_{+}},n{{W}_{-}} \right\}}{{{Z}^{m}}({{W}_{+}}){{Z}^{n}}({{W}_{-}})}.\end{eqnarray} \tag{ C.3 }$$

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

$$\begin{eqnarray}&&Z(W)=(\pm )\sqrt{{\rm det} \left( 2{\rm cosh} \frac{W}{2} \right)},\qquad \left\{ {{W}_{1}},{{W}_{2}} \right\}=(\pm )\sqrt{{\rm det} (1+{{{\rm e}}^{{{W}_{1}}}}{{{\rm e}}^{{{W}_{2}}}})},\end{eqnarray} \tag{ 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 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 $O_{+}^{3}$ and $O_{+}^{2}{{O}_{-}}$, 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 ${{\Gamma }_{\pm }}$ with the result

$$\begin{eqnarray}&&{\rm Tr}(O_{+}^{3})=\pm \sqrt{{\rm det} \left( \frac{1+3\Gamma _{+}^{2}}{4} \right)},\qquad {\rm Tr}(O_{+}^{2}{{O}_{-}})=\sqrt{{\rm det} \left( \frac{1+\Gamma _{+}^{2}+2{{\Gamma }_{+}}{{\Gamma }_{-}}}{4} \right)}.\end{eqnarray} \tag{ C.5 }$$

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 ${{\Gamma }_{+}}$ contains a quadruplet of purely imaginary eigenvalues. For the Ising chain, this can happen only in the symmetry-broken phase, $h\lt 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.