Twisted and untwisted negativity spectrum of free fermions

TL;DR

This paper investigates the negativity spectrum of free fermions under different definitions of partial transpose, revealing complex eigenvalue distributions and extending entanglement diagnostics to fermionic systems.

Contribution

It introduces and analyzes two definitions of fermionic partial transpose, deriving their spectral properties and eigenvalue distributions using path-integral formalism.

Findings

Two different fermionic PT definitions lead to Hermitian and pseudo-Hermitian density matrices.

Analytical expressions for the moments of the spectra are derived.

Numerical results confirm the analytical predictions.

Abstract

A basic diagnostic of entanglement in mixed quantum states is known as the positive partial transpose (PT) criterion. Such criterion is based on the observation that the spectrum of the partially transposed density matrix of an entangled state contains negative eigenvalues, in turn, used to define an entanglement measure called the logarithmic negativity. Despite the great success of logarithmic negativity in characterizing bosonic many-body systems, generalizing the operation of PT to fermionic systems remained a technical challenge until recently when a more natural definition of PT for fermions that accounts for the Fermi statistics has been put forward. In this paper, we study the many-body spectrum of the reduced density matrix of two adjacent intervals for one-dimensional free fermions after applying the fermionic PT. We show that in general there is a freedom in the definition of such operation which leads to two different definitions of PT: the resulting density matrix is Hermitian in one case, while it becomes pseudo-Hermitian in the other case. Using the path-integral formalism, we analytically compute the leading order term of the moments in both cases and derive the distribution of the corresponding eigenvalues over the complex plane. We further verify our analytical findings by checking them against numerical lattice calculations.

Equations

$$S({\rho }_A)=-\text{Tr}({\rho }_A\ln {\rho }_A)=-\sum _i{\alpha }_i^2\ln {\alpha }_i^2,$$

$${\mathcal{R}}_n({\rho }_A)=\frac{1}{1-n}\ln \text{Tr}({\rho }_A^n)=\frac{1}{1-n}\ln \sum _i{\alpha }_i^{2n},$$

$${\rho }^{T_A}=\sum _i{\rho }_{ijkl}|e_A^{(k)},e_B^{(j)}⟩⟨e_A^{(i)},e_B^{(l)}|.$$

$$\mathcal{N}(\rho )$$

$$\mathcal{E}(\rho )$$

$${\mathcal{N}}_n(\rho )=\ln \text{Tr}({\rho }^{T_A}{)}^n,$$

$$\mathcal{E}(\rho )=\underset{n\to 1/2}{\lim }{\mathcal{N}}_{2n}(\rho ).$$

$$P(\lambda )=\sum _i\delta (\lambda -{\lambda }_i)$$

$$\hat{H}=-\sum _j[t(f_{j+1}^{†}{f}_j+\text{H.c.})+\mu f_j^{†}{f}_j],$$

$$c_{2j-1}:=f_j^{†}+f_j,c_{2j}:=i(f_j-f_j^{†}),j=1,\dots ,N.$$

$$X=\sum _{k=1}^{2N}\sum _{p_1<p_2\cdots <p_k}{X}_{p_1\cdots p_k}{c}_{p_1}\cdots c_{p_k},$$

$$\rho =\sum _{k_1,k_2}^{k_1+k_2=\mathrm{even}}{\rho }_{p_1\cdots p_{k_1},q_1\cdots q_{k_2}}{a}_{p_1}\cdots a_{p_{k_1}}{b}_{q_1}\cdots b_{q_{k_2}},$$

$${\rho }^{T_A}:=\sum _{k_1,k_2}^{k_1+k_2=\mathrm{even}}{\rho }_{p_1\cdots p_{k_1},q_1\cdots q_{k_2}}{i}^{k_1}{a}_{p_1}\cdots a_{p_{k_1}}{b}_{q_1}\cdots b_{q_{k_2}},$$

$$({\rho }^{T_A}{)}^{†}$$

$$({\rho }^{T_A}{)}^{T_A}$$

$${\rho }^{{\tilde{T}}_A}:={\rho }^{T_A}(-1{)}^{F_A}.$$

$$({\rho }^{{\tilde{T}}_A}{)}^{{\tilde{T}}_A}=\rho ,$$

$${\mathcal{N}}_n^{(\mathrm{ns})}(\rho )=\ln \text{Tr}[({\rho }^{T_A}{)}^n],{\mathcal{N}}_n^{(\mathrm{r})}(\rho )=\ln \text{Tr}[({\rho }^{{\tilde{T}}_A}{)}^n],$$

$$\mathcal{E}(\rho )=\underset{n\to 1/2}{\lim }{\mathcal{N}}_{2n}^{(\mathrm{r})}.$$

$$R_n:=\text{Tr}[{\mathcal{O}}^n].$$

$$f(s):=\frac{1}{\pi }\sum _{n=1}^{\infty }{R}_n{s}^{-n}=\frac{1}{\pi }\int d\lambda \frac{\lambda P(\lambda )}{s-\lambda }.$$

$$P(\lambda )=\frac{1}{\lambda }\underset{ϵ\to 0}{\lim }\text{Im}f(\lambda -iϵ).$$

$$n(\lambda )={\int }_{\lambda }^{{\lambda }_{\text{max}}}d\lambda P(\lambda ).$$

$$R_n=r_n\exp \left(-bn+\frac{a}{n}\right),\forall n,$$

$$P(\lambda ;a,b)$$

$$n(\lambda ;a,b)$$

$$\text{Spec}[{\rho }^{T_A}]$$

$$\text{Spec}[{\rho }^{{\tilde{T}}_A}]$$

$${\rho }_{\Omega }=\frac{1}{\mathcal{Z}(\Omega )}\exp \left(\frac{1}{4}\sum _{j,k=1}^{2N}{\Omega }_{jk}{c}_j{c}_k\right),$$

$${\Gamma }_{jk}=\frac{1}{2}\text{Tr}({\rho }_{\Omega }[c_j,c_k]),$$

$$\Gamma =\tanh \left(\frac{\Omega }{2}\right),e^{\Omega }=\frac{\mathbb{I}+\Gamma }{\mathbb{I}-\Gamma }.$$

$${\mathcal{R}}_n(\rho )=\frac{1}{1-n}\sum _{j=1}^N\ln \left[{\left(\frac{1-{\nu }_j}{2}\right)}^n+{\left(\frac{1+{\nu }_j}{2}\right)}^n\right].$$

Full text

Twisted and untwisted negativity spectrum of free fermions

Hassan Shapourian1, Paola Ruggiero2, Shinsei Ryu1, Pasquale Calabrese2,3

1 James Franck Institute and Kadanoff Center for Theoretical Physics, University of Chicago, IL 60637

2 SISSA and INFN, via Bonomea 265, 34136 Trieste, Italy

3 International Centre for Theoretical Physics (ICTP), I-34151, Trieste, Italy

Abstract

A basic diagnostic of entanglement in mixed quantum states is known as the positive partial transpose (PT) criterion. Such criterion is based on the observation that the spectrum of the partially transposed density matrix of an entangled state contains negative eigenvalues, in turn, used to define an entanglement measure called the logarithmic negativity. Despite the great success of logarithmic negativity in characterizing bosonic many-body systems, generalizing the operation of PT to fermionic systems remained a technical challenge until recently when a more natural definition of PT for fermions that accounts for the Fermi statistics has been put forward. In this paper, we study the many-body spectrum of the reduced density matrix of two adjacent intervals for one-dimensional free fermions after applying the fermionic PT. We show that in general there is a freedom in the definition of such operation which leads to two different definitions of PT: the resulting density matrix is Hermitian in one case, while it becomes pseudo-Hermitian in the other case. Using the path-integral formalism, we analytically compute the leading order term of the moments in both cases and derive the distribution of the corresponding eigenvalues over the complex plane. We further verify our analytical findings by checking them against numerical lattice calculations.

Contents

1. 1 Introduction

2. 2 Preliminary remarks

3. 2.1 Twisted and untwisted partial transpose for fermions

4. 2.2 The moment problem

5. 2.3 Partial transpose of Gaussian states

6. 3 Spacetime picture for the moments of partial transpose

7. 4 The spectrum of partial transpose

8. 4.1 Spectrum of ${\rho^{T_{A}}}$

9. 4.1.1 Adjacent intervals

10. 4.1.2 Bipartite geometry

11. 4.2 Spectrum of $\rho^{\widetilde{T}_{A}}$

12. 4.2.1 Adjacent intervals

13. 4.2.2 Bipartite geometry

14. 5 Conclusions

15. A Twist fields, bosonization, etc.

16. B Rényi negativity for disjoint intervals

17. B.1 Moments of $\rho^{{T}_{A}}$

18. B.2 Moments of $\rho^{\widetilde{T}_{A}}$

19. C Partial transpose with supersymmetric trace

20. D Negativity of bosonic scalar field theory

1 Introduction

Entanglement is an intrinsic property of quantum systems beyond classical physics. Having efficient frameworks to compute entanglement between two parts of a system is essential not only for fundamental interests such as characterizing phases of matter [1, 2, 3, 4] and spacetime physics [5] but also for application purposes such as identifying useful resources to implement quantum computing processes. For a bipartite Hilbert space $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, it is easy to determine whether a pure state $\ket{\Psi}$ is entangled or not: a product state, i.e. any state of the form $\ket{\Phi_{A}}\otimes\ket{\Phi_{B}}$, is unentangled (separable), while a superposition state $\ket{\Psi}=\sum_{i}\alpha_{i}\ket{\Phi_{A}^{(i)}}\otimes\ket{\Phi_{B}^{(i)}}$, where $\ket{\Phi_{A/B}^{(i)}}$ is a set of local orthogonal states, is entangled. The amount of entanglement in a given state can be quantified by the entropy of information within either subsystem $A$ or $B$, in the form of the von Neumann entropy

[TABLE]

or the Rényi entanglement entropies (REEs)

[TABLE]

where $\rho_{A}=\text{Tr}_{B}(\ket{\Psi}\bra{\Psi})=\sum_{i}\alpha_{i}^{2}\ket{\Phi_{A}^{(i)}}\bra{\Phi_{A}^{(i)}}$ is the reduced density matrix acting on $\mathcal{H}_{A}$. Notice that $S(\rho_{A})=S(\rho_{B})$ and ${\cal R}_{n}(\rho_{A})={\cal R}_{n}(\rho_{B})$ and clearly, $S,{\cal R}_{n}\geq 0$ where the equality holds for a product state. For analytical calculations, $S$ is usually obtained from ${\cal R}_{n}$ via $S=\lim_{n\to 1}{\cal R}_{n}$.

It is well-known that eigenvalues of density matrices, i.e. the entanglement spectrum, contains more information than merely the entanglement entropies. The entanglement spectrum has been studied and utilized toward better understanding of the phases of matter [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23], broken-symmetry phases [24, 25, 26, 27, 28], and more exotic phases such as many-body localized states [29, 30, 31]. In particular, in the context of conformal field theories (CFTs) in (1+1)d the distribution of eigenvalues was analytically derived [32] and was shown to obey a universal scaling function which depends only on the central charge of the underlying CFT. The obtained scaling function for the distribution of the entanglement spectrum at criticality was further substantiated numerically [33, 34], especially for matrix product state representation at critical points [35, 36, 37].

It turned out that extending the above ideas to mixed states where the system is described by a density matrix $\rho$ is not as easy as it may seem. A product state $\rho=\rho_{A}\otimes\rho_{B}$ is similarly unentangled. However, a large class of states, called separable states, in the form of $\rho=\sum_{i}p_{i}\rho_{A}^{(i)}\otimes\rho_{B}^{(i)}$ with $p_{i}\geq 0$ are classically correlated and do not contain any amount of entanglement. Hence, the fact that superposition implies entanglement in pure states does not simply generalize to the entanglement in mixed states. The positive partial transpose (PPT) [38, 39, 40, 41, 42, 43, 44] is a test designed to diagnose separable states based on the fact that density matrices are positive semi-definite operators. The PT of a density matrix $\rho=\sum_{i}\rho_{ijkl}\ket{e_{A}^{(i)},e_{B}^{(j)}}\bra{e_{A}^{(k)},e_{B}^{(l)}}$ written in a local orthonormal basis $\{\ket{e_{A}^{(k)}},\ket{e_{B}^{(j)}}\}$ is defined by exchanging the indices of subsystem $A$ (or $B$) as in

[TABLE]

Note that $\rho^{T_{A}}$ is a Hermitian operator and the PPT test follows by checking whether or not $\rho^{T_{A}}$ contains any negative eigenvalue. A separable state passes the PPT test, i.e. all the eigenvalues of $\rho^{T_{A}}$ are non-negative, whereas an inseparable (i.e., entangled) state yields negative eigenvalues after PT111A technical point is that there exists a set of inseparable states which also pass the PPT test [45]. They are said to contain bound entanglement which cannot be used for quantum computing processes such as teleportation [46]. This issue is beyond the scope of our paper and we do not elaborate further here.. Hence, the PPT criterion can be used to decide whether a given density matrix is separable or not. Similar to the entropic measures of pure-state entanglement in (1) and (2), the (logarithmic) entanglement negativity associated with the spectrum of the partially transposed density matrix is defined as a candidate to quantify mixed-state entanglement [47, 48, 49],

[TABLE]

where $\left\lVert A\right\rVert=\text{Tr}\sqrt{AA^{\dagger}}$ is the trace norm. When $A$ is Hermitian, the trace norm is simplified into the sum of the absolute value of the eigenvalues of $A$. Hence, the above quantities measure the negativity of the eigenvalues of $\rho^{T_{A}}$. It is also useful to define the moments of the PT (aka Rényi negativity) via

[TABLE]

where the logarithmic negativity is obtained from analytic continuation

[TABLE]

Note that for a pure state $\rho=\ket{\Psi}\bra{\Psi}$, we have ${\cal E}(\rho)={\cal R}_{1/2}(\rho_{A})$ where $\rho_{A}$ is the reduced density matrix on $\mathcal{H}_{A}$. The entanglement negativity has been used to characterize mixed states in various quantum systems such as in harmonic oscillator chains [50, 51, 52, 53, 54, 55, 56, 57, 58], quantum spin models [59, 60, 61, 62, 63, 64, 65, 66, 67, 68], (1+1)d conformal and integrable field theories [69, 70, 71, 72, 73, 74, 75, 76], topologically ordered phases of matter in (2+1)d [77, 78, 79, 80, 81], and in out-of-equilibrium situations [82, 83, 84, 85, 86, 87], as well as holographic theories [88, 89, 90, 91, 92, 93] and variational states [94, 95, 96, 97]. Moreover, the PT was used to construct topological invariants for symmetry protected topological (SPT) phases protected by anti-unitary symmetries [98, 99, 100, 101] and there are experimental proposals to measure it with cold atoms [102, 103].

Unlike the entanglement spectrum which has been studied extensively, less is known about the spectrum of partially transposed density matrices in many-body systems. It is true that the PPT test which predates the entanglement spectrum is based on the eigenvalues of the PT, but the test only uses the sign of the eigenvalues. Therefore, studying the spectrum of the PT could be useful in characterizing quantum phases of matter. Furthermore, the fact that PT is applicable to extract entanglement at finite-temperature states and that the eigenvalues have a sign structure (positive/negative) may help unravel some new features beyond the entanglement spectrum. Recently, the distribution of eigenvalues of the PT, dubbed as the negativity spectrum, was studied for CFTs in (1+1)d [72]. It was found that the negativity spectrum is universal and depends only on the central charge of the CFT, similar to the entanglement spectrum, while the precise form of the spectrum depends on the sign of the eigenvalues. This dependence is weak for bulk eigenvalues, whereas it is strong at the spectrum edges.

In this paper, we would like to study the negativity spectrum in fermionic systems. The PT of fermionic density matrices however involves some subtleties due to the Fermi statistics (i.e., anti-commutation relation of fermion operators). Initially, a procedure for the PT of fermions based on the fermion-boson mapping (Jordan-Wigner transformation) was proposed [104] and was also used in the subsequent studies [105, 106, 107, 108, 109, 110, 111]. However, this definition turned out to cause certain inconsistencies within fermionic theories such as violating the additivity property and missing some entanglement in topological superconductors, and give rise to incorrect classification of time-reversal symmetric topological insulators and superconductors. Additionally, according to this definition it is computationally hard to find the PT (and calculate the entanglement negativity) even for free fermions, since the PT of a fermionic Gaussian state is not Gaussian. This motivates us to use another way of implementing a fermionic PT which was proposed recently by some of us in the context of time-reversal symmetric SPT phases of fermions [100, 101, 112]. This definition does not suffer from the above issues and at the same time the associated entanglement quantity is an entanglement monotone [113]. From a practical standpoint, the latter definition has the merit that the partially transposed Gaussian state remains Gaussian and hence can be computed efficiently for free fermions. A detailed survey of differences between the two definitions of PT from both perspectives of quantum information and condensed matter theory (specifically, topological phases of fermions) is discussed in Refs. [112, 113].

Before we get into details of the fermionic PT in the coming sections, let us finish this part with a summary of our main findings. We study the distribution of the many-body eigenvalues $\lambda_{i}$ of the partially transposed reduced density matrix,

[TABLE]

for one-dimensional free fermions. As a lattice realization, we consider the hopping Hamiltonian on a chain

[TABLE]

where the fermion operators $f_{j}$ and $f_{j}^{\dagger}$ obey the anti-commutation relation $\{f_{i},{f_{j}}^{\dagger}\}=f_{i}{f_{j}}^{\dagger}+{f_{j}}^{\dagger}f_{i}=\delta_{ij}$ and $\{f_{i},f_{j}\}=\{{f_{i}}^{\dagger},{f_{j}}^{\dagger}\}=0$.

Recall that using the regular (matrix) PT – we will refer to it as the bosonic PT –, which applies to generic systems where local operators commute, the obtained PT density matrix is a Hermitian operator and its eigenvalues are either negative or positive. However, it turned out that for fermions a consistent definition of PT involves a phase factor as we exchange indices in (3) and in general one can define two types of PT operation. As we will explain in detail, these two types correspond to the freedom of spacetime boundary condition for fermions associated with the fermion-number parity symmetry. We reserve ${\rho^{T_{A}}}$ and $\rho^{\widetilde{T}_{A}}$ to denote the fermionic PT which leads to anti-periodic (untwisted) and periodic (twisted) boundary conditions along fundamental cycles of the spacetime manifold, respectively. We should note that ${\rho^{T_{A}}}$ is pseudo-Hermitian222A pseudo-Hermitian operator $H$ is defined by $\eta H^{\dagger}\eta^{-1}=H$ with $\eta^{2}=1$ where $\eta$ is a unitary Hermitian operator satisfying $\eta^{\dagger}\eta=\eta\eta^{\dagger}=1$ and $\eta=\eta^{\dagger}$. Essentially, pseudo-Hermiticity is a generalization of Hermiticity, in that it implies Hermiticity when $\eta=1$. and may contain complex eigenvalues, while $\rho^{\widetilde{T}_{A}}$ is Hermitian and its eigenvalues are real. We use the spacetime path integral formulation to analytically calculate the negativity spectrum. In the case of $\rho^{\widetilde{T}_{A}}$, we obtain results very similar to those of previous CFT work [72], where the distribution of positive and negative eigenvalues are described by two universal functions. In the case of ${\rho^{T_{A}}}$, we observe that the eigenvalues are complex but they have a pattern and fall on six branches in complex plane with a quantized complex phase of $\angle\lambda=2\pi n/6$. We show that the spectrum is reflection symmetric with respect to the real axis and the eigenvalue distributions are described by four universal functions along $\angle\lambda=0,\pm 2\pi/6,\pm 4\pi/6,\pi$ branches. We further verify our findings by checking them against numerical lattice simulations.

The rest of our paper is organized as follows: in Section 2 we provide a brief review of partial transpose for fermions, in Section 3 we discuss the spacetime path-integral formulation of the moments of partially transposed density matrices. The spectrum of the twisted and untwisted partial transpose is analytically derived in Section 4 for different geometries, where numerical checks with free fermions on the lattice are also provided. We close our discussion by some concluding remarks in Section 5. In several appendices, we give further details of the analytical calculations and make connections with other related concepts.

2 Preliminary remarks

In this section, we review some basic materials which we use in the next sections: the definition of PT for fermions, how to extract the distribution of the eigenvalues of an operator from its moments, and some properties of partially transposed Gaussian states.

2.1 Twisted and untwisted partial transpose for fermions

In this part, we briefly discuss some background materials on our definitions of PT for fermions. More details can be found in Refs. [112, 113]. We consider a fermionic Fock space ${\cal H}$ generated by $N$ local fermionic modes $f_{j}$, $j=1,\cdots,N$. The Hilbert space is spanned by $\ket{n_{1},n_{2},\cdots,n_{N}}$ which is a string of occupation numbers $n_{j}=0,1$. We define the Majorana (real) fermion operators in terms of canonical operators as

[TABLE]

These operators satisfy the commutation relation $\{c_{j},c_{k}\}=2\delta_{jk}$ and generate a Clifford algebra. Any operator $X$ acting on $\mathcal{H}$ can be expressed in terms of a polynomial of $c_{j}$’s,

[TABLE]

where $X_{p_{1}\dots p_{k}}$ are complex numbers and fully antisymmetric under permutations of $\{1,\dots,k\}$. A density matrix has an extra constraint, i.e., it commutes with the total fermion-number parity operator, $[\rho,(-1)^{F}]=0$ where $F=\sum_{j}f_{j}^{\dagger}f_{j}$. This constraint entails that the Majorana operator expansion of $\rho$ only contains even number of Majorana operators, i.e., $k$ in the above expression is even.

To study the entanglement, we consider a bipartite Hilbert space $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ spanned by $f_{j}$ with $j=1,\cdots,N_{A}$ in subsystem $A$ and $j=N_{A}+1,\cdots,N_{A}+N_{B}$ in subsystem $B$. Then, a generic density matrix on $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ can be expanded in Majorana operators as

[TABLE]

where $\{a_{j}\}$ and $\{b_{j}\}$ are Majorana operators acting on $\mathcal{H}^{A}$ and $\mathcal{H}^{B}$, respectively, and the even fermion-number parity condition is indicated by the condition $k_{1}+k_{2}={\rm even}$. Our definition of the PT for fermions is given by [112, 101]

[TABLE]

and similarly for $\rho^{T_{B}}$. It is easy to see that the subsequent application of the PT with respect to the two subsystems leads to the full transpose $(\rho^{T_{A}})^{T_{B}}=\rho^{T}$, i.e. reversing the order of Majorana fermion operators. In addition, the definition (13) implies that

[TABLE]

where $(-1)^{F_{A}}$ is the fermion-number parity operator on $H_{A}$, i.e. $F_{A}=\sum_{j\in A}f_{j}^{\dagger}f_{j}$. The first identity, namely the pseudo-Hermiticity, can be understood as a consequence of the fact that $({\rho^{T_{A}}})^{\dagger}$ is defined the same as (13) by replacing $i^{k_{1}}$ with $(-i)^{k_{1}}$. The second identity reflects the fact that the fermionic PT is related to the action of time-reversal operator of spinless fermions in the Euclidean spacetime [101]. We should note that the matrix resulting from the PT is not necessarily Hermitian and may have complex eigenvalues, although $\text{Tr}{\rho^{T_{A}}}=1$. The existence of complex eigenvalues is a crucial property which was used in the context of SPT invariants to show that the complex phase of $\text{Tr}(\rho{\rho^{T_{A}}})$, which represents a partition function on a non-orientable spacetime manifold, is a topological invariant. For instance, $\text{Tr}(\rho{\rho^{T_{A}}})=e^{i2\pi\nu/8}$ for time-reversal symmetric topological superconductors (class BDI) which implies the $\mathbb{Z}_{8}$ classification. (Here $\nu\in\mathbb{Z}_{8}$ is the topological invariant). Nevertheless, we may still use Eq. (5) to define an analog of entanglement negativity for fermions and calculate the trace norm in terms of square root of the eigenvalues of the composite operator $\rho_{\times}=[({\rho^{T_{A}}})^{\dagger}{\rho^{T_{A}}}]$, which is a Hermitian operator with real positive eigenvalues. On the other hand, from Eq. (14) we realize that $\rho_{\times}=(\rho^{\widetilde{T}_{A}})^{2}$ where we introduce the twisted PT by

[TABLE]

It is easy to see from Eq. (14) that this operator is Hermitian and then similar to the bosonic PT always contains real eigenvalues. It is worth noting that

[TABLE]

in contrast with the untwisted PT (15). As we will see shortly, this difference between ${\rho^{T_{A}}}$ and $\rho^{\widetilde{T}_{A}}$ in the operator formalism will show up as anti-periodic and periodic boundary conditions across the fundamental cycles of spacetime manifold in the path-integral formalism. The central result of our paper is to report analytical results for the spectrum of ${\rho^{T_{A}}}$ and $\rho^{\widetilde{T}_{A}}$.

2.2 The moment problem

In the replica approach to logarithmic negativity (5) and negativity spectrum, one first has to calculate the moments of PT, aka Rényi negativity (RN),

[TABLE]

which are fermionic counterparts of the bosonic definition in Eq. (6). The superscripts $(\textrm{ns})$ and $(\textrm{r})$ stand for Neveu-Schwarz and Ramond respectively (the reason for this will be clear from the path integral representation of such quantities, see Section 3 below). Thus, the analog of analytic continuation (7) to obtain the logarithmic negativity is

[TABLE]

In the following, we review a general framework to analytically obtain the distribution of eigenvalues of density matrix (or its transpose) from the moments. This method was originally used to derive the entanglement spectrum of (1+1)d CFTs [32]. Suppose we have an operator ${\cal O}$ whose moments are of the form

[TABLE]

In terms of the eigenvalues of ${\cal O}$, $\{\lambda_{j}\}$, we have $R_{n}=\sum_{j}\lambda_{j}^{n}=\int P(\lambda)\lambda^{n}d\lambda$, where $P(\lambda)$ is the associated distribution function (see Eq. (8)). The goal is to find $P(\lambda)$ by making use of the specific form of $R_{n}$ in (20). The essential idea is to compute the Stjilties transform

[TABLE]

Assuming that the eigenvalues are real, the distribution function can be easily read off from the relation

[TABLE]

In the following we are going to focus on the complementary cumulative distribution function or simply the tail distribution, being a very simple object to be accessed for numerical comparison

[TABLE]

For specific types of operators such as the density matrices and their PT in (1+1)d CFTs, the moments can be cast in the form,

[TABLE]

where $a,b\in\mathbb{R}$, $b>0$ and $r_{n}$ are non-universal constant. In such cases, the distribution function is found to be [63]

[TABLE]

and the corresponding tail distribution is given by

[TABLE]

where $J_{\alpha}(x)$ and $I_{\alpha}(x)$ are the regular Bessel functions and modified Bessel functions of the first kind, respectively. Note that (27) and (30) are derived by ignoring the presence of the constants $r_{n}$ in (24). This relies on the assumption that they do not change significantly upon varying $n$, i.e., $\lim_{n\to\infty}\frac{1}{n}\ln r_{n}<\infty$. The very same assumption has been adopted for the entanglement and bosonic negativity spectrum in Refs. [32] and [72] where the derived distribution functions agree with the numerically obtained spectra.

2.3 Partial transpose of Gaussian states

Here, we discuss how to compute the spectrum of the PT of a Gaussian state from the corresponding covariance matrix. The idea is similar to that of the entanglement spectrum, while there are some differences as the covariance matrix associated with the partially transposed density matrix may contain complex eigenvalues. Before we continue, let us summarize the structure of the many-body spectrum of ${\rho^{T_{A}}}$ and $\rho^{\widetilde{T}_{A}}$ for free fermions,

[TABLE]

where repeating values mean two fold degeneracy. We should note that the pseudo-Hermiticity of ${\rho^{T_{A}}}$ (14) ensures that the complex-valued subset of many-body eigenvalues of ${\rho^{T_{A}}}$ appear in complex conjugate pairs. This property is general and applicable to any density matrix beyond free fermions. An immediate consequence of this property is that any moment of ${\rho^{T_{A}}}$ is guaranteed to be real-valued.

A Gaussian density matrix in the Majorana fermion basis (10) is defined by

[TABLE]

where $\Omega$ is a pure imaginary antisymmetric matrix and ${\cal Z}(\Omega)=\pm\sqrt{\det\left(2\cosh\frac{\Omega}{2}\right)}$ is the normalization constant. We should note that the spectrum of $\Omega$ is in the form of $\pm\omega_{j},\ j=1,\dots,N$ and the $\pm$ sign ambiguity in ${\cal Z}(\Omega)$ is related to the square root of determinant where we need to choose one eigenvalue for every pair $\pm\omega_{j}$. The sign is fixed by the Pfaffian. This density matrix can be uniquely characterized by its covariance matrix,

[TABLE]

which is a $2N\times 2N$ matrix. These two matrices are related by

[TABLE]

Furthermore, one can consider a generic Gaussian operator which is also defined through Eq. (38), but without requiring that the spectrum is pure imaginary. An equivalent description in terms of the covariance matrix is also applicable for such operators. The only difference is that the eigenvalues do not need to be real. Let us recall how Rényi entropies (2) are computed for Gaussian states. The density matrix (38) can be brought into a diagonal form $\rho_{\Omega}={\cal Z}^{-1}\exp\left(\frac{i}{2}\sum_{n}\omega_{n}d_{2n}d_{2n-1}\right)$, where $\omega_{n}$ is obtained from an orthogonal transformation of $\Omega$. In terms of the eigenvalues of $\Gamma$, denoted by $\pm\nu_{j}$, we have $\rho_{\Omega}=\prod_{n}(1+i\nu_{n}d_{2n}d_{2n-1})/2$, leading to

[TABLE]

We consider a density matrix on a bipartite Hilbert space (12) where the covariance matrix takes a block matrix form as

[TABLE]

Here, $\Gamma_{AA}$ and $\Gamma_{BB}$ denote the reduced covariance matrices of subsystems $A$ and $B$, respectively; while $\Gamma_{AB}=\Gamma_{BA}^{\dagger}$ describes the correlations between them. We define the covariance matrix associated with a partially transposed Gaussian state by

[TABLE]

where $[\Gamma_{+}]_{ij}=\frac{1}{2}\text{Tr}({\rho^{T_{A}}}[c_{i},c_{j}])$ and $[\Gamma_{-}]_{ij}=\frac{1}{2}\text{Tr}({\rho^{T_{A}}}^{\dagger}[c_{i},c_{j}])$. We should note that $\Gamma_{+}$ and $\Gamma_{-}$ have identical eigenvalues while they do not necessarily commute $[\Gamma_{+},\Gamma_{-}]\neq 0$. In general, the eigenvalues of $\Gamma_{+}$ appear in quartets $(\pm\nu_{k},\pm\nu_{k}^{\ast})$ when $\text{Re}[\nu_{k}]\neq 0$ and $\text{Im}[\nu_{k}]\neq 0$ or doublet $\pm\nu_{k}$ when $\text{Re}[\nu_{k}]=0$ (i.e., pure imaginary) or $\text{Im}[\nu_{k}]=0$ (i.e., real) . $\pm$ is because of skew symmetry $\Gamma^{T}_{\pm}=-\Gamma_{\pm}$. In addition, the pseudo-Hermiticity of ${\rho^{T_{A}}}$ (14) implies that

[TABLE]

where $U_{1}=(-\mathbb{I}_{A}\oplus\mathbb{I}_{B})$ is the matrix associated with the operator $(-1)^{F_{A}}$. This means that for every eigenvalue $\nu_{k}$ its complex conjugate $\nu_{k}^{\ast}$ is also an eigenvalue. As a result, the moments of PT can be written as

[TABLE]

Note that the sum is now over half of the eigenvalues (say in the upper half complex plane), due to the structure discussed above.

For $\rho^{\widetilde{T}_{A}}$ we use the multiplication rule for the Gaussian operators where the resulting Gaussian matrix is given by

[TABLE]

which is manifestly Hermitian due to the identity (48). Using Eq. (38), the normalization factor is found to be ${\cal Z}_{\widetilde{T}_{A}}=\text{Tr}(\rho^{\widetilde{T}_{A}})=\text{Tr}[\rho(-1)^{F_{A}}]=\sqrt{\det\Gamma_{AA}}$. From (40) we construct the covariance matrix $\widetilde{\Gamma}_{\pm}=\tanh(\widetilde{\Omega}/2)$ and compute the moments of $\rho^{\widetilde{T}_{A}}$ by

[TABLE]

where $\pm\tilde{\nu}_{j}$ are eigenvalues of $\widetilde{\Gamma}_{\pm}$ which are guaranteed to be real. Consequently, the logarithmic negativity (19) is given by

[TABLE]

For particle-number conserving systems such as the lattice model in (9), the covariance matrix is simplified into the form $\Gamma=\sigma_{2}\otimes\gamma$ where $\gamma=(\mathbb{I}-2C)$ and $C_{ij}=\text{Tr}(\rho f_{i}^{\dagger}f_{j})$ is the correlation matrix and $\sigma_{2}$ is the second Pauli matrix acting on the even/odd indices of Majorana operators $(c_{2j},c_{2j-1})$. In this case, the transformed correlation matrix for ${\rho^{T_{A}}}$ is given by

[TABLE]

The eigenvalues can be divided to two categories: complex eigenvalues $\nu_{k}$, $\text{Im}[\nu_{k}]\neq 0$ and real eigenvalues $u_{k}$, $\text{Im}[u_{k}]=0$. The pseudo-Hermiticity property leads to the identity $\gamma_{\pm}^{\dagger}=U_{1}\gamma_{\pm}U_{1}$ which implies that complex eigenvalues appear in pairs $(\nu_{k},\nu_{k}^{\ast})$. Therefore, the many-body eigenvalues follow the form,

[TABLE]

where $\bm{\sigma}=\{\sigma_{k}=\pm\}$ is a string of signs. Clearly, the many-body eigenvalues appear in two categories as well: complex conjugate pairs $(\lambda_{j},\lambda_{j}^{\ast})$ and real eigenvalues which are not necessarily degenerate.

We can also derive a simple expression for the correlation matrix $\widetilde{C}=(\mathbb{I}-\widetilde{\gamma})/2$ associated with $\rho^{\widetilde{T}_{A}}$,

[TABLE]

3 Spacetime picture for the moments of partial transpose

In the following two sections, we compute the moments of the partially transposed density matrix and ultimately the logarithmic negativity. First, we develop a general method using the replica approach [69, 70, 110] and provide an equivalent spacetime picture of the Rényi negativity.

Before we proceed, let us briefly review the replica approach to find the entanglement entropy. Next, we make connections to our construction of PT. A generic density matrix can be represented in the fermionic coherent state as

[TABLE]

where $\alpha$, $\bar{\alpha}$, $\beta$ and $\bar{\beta}$ are independent Grassmann variables and we omit the real-space (and possibly other) indices for simplicity. The trace formula then reads

[TABLE]

where the subscripts in $\psi_{i}$ and $\bar{\psi}_{i}$ denote the replica indices and ${T}$ is called the twist matrix,

[TABLE]

The above expression can be viewed as a partition function on a $n$-sheet spacetime manifold where the $n$ flavors (replicas) $\psi_{i}$ are glued in order along the cuts. Alternatively, one can consider a multi-component field $\Psi=(\psi_{1},\cdots,\psi_{n})^{T}$ on a single-sheet spacetime. This way when we traverse a close path through the interval the field gets transformed as $\Psi\mapsto T\Psi$. Hence, each interval can be represented by two branch points ${\cal T}_{n}$ and ${\cal T}_{n}^{-1}$ –the so-called twist fields– and the REE of one interval can be written as a two-point correlator [114],

[TABLE]

where $u$ and $v$ denote the real space coordinates of the two ends of the interval defining the subsystem $A$.

Let us now derive analogous relations for the moments of partially transposed density matrix. Using the definition of the PT in the coherent state basis [112]

[TABLE]

we write the general expression for the moments of ${\rho^{T_{A}}}$ as

[TABLE]

where $\psi_{js}$ and $\bar{\psi}_{js}$ refer to the field defined within the $s=A,B$ interval of $j$th replica. Here, we are dealing with two intervals where the twist matrices are $T$ and $T^{-1}$ as shown in Fig. 1(a). Therefore, it can be written as a four-point correlator (1(b))

[TABLE]

Note that the order of twist fields are reversed for the first interval.

From the coherent state representation, we can also write the moments of $\rho^{\widetilde{T}_{A}}$

[TABLE]

The twist matrix for interval A is modified to be

[TABLE]

which can be viewed as a gauge transformed twist matrix $T^{-1}$. Analogously, Eq. (71) can be written in terms of a four-point correlator

[TABLE]

where $\widetilde{\cal T}_{n}$ and $\widetilde{\cal T}^{-1}_{n}$ are twist fields associated with $\widetilde{T}$.

For fermions with a global $U(1)$ gauge symmetry (i.e., particle-number conserving systems) there is a freedom to twist boundary condition along the fundamental cycles (e.g. the dashed-line path in Fig. 1(a)) of the spacetime manifold by a $U(1)$ phase (or holonomy). The boundary conditions are independent and in principle can be different for different pairs of sheets. If we assume a replica symmetry (i.e. uniform boundary conditions) $\psi_{i}\mapsto e^{i\alpha}\psi_{i}$, the expression for the PT moments in the operator formalism is given by

[TABLE]

Let us mention that some related quantities such as $\text{Tr}[(\rho\ e^{i\alpha F})^{n}]$ were previously introduced and dubbed charged entanglement entropies [115]. They were further used to determine symmetry resolved entanglement entropies which is the contribution from the density matrix to the entanglement entropies when projected onto a given particle-number sector [116, 117].

From (78), we get a family of RN parametrized by $\alpha$. However, for a generic fermionic system (including superconductors), the $U(1)$ symmetry is reduced to $\mathbb{Z}_{2}$ fermion-parity symmetry. Hence, the two quantities of general interest would be

[TABLE]

We should reemphasize that either quantities are described by a partition function on the same spacetime manifold (Fig. 1) as in the case of bosonic systems [70], while they differ in the boundary conditions for fundamental cycles of the manifold. In other words, ${Z}_{{\cal N}_{n}}^{({\rm ns})}$ and ${Z}_{{\cal N}_{n}}^{({\rm r})}$ correspond to anti-periodic (i.e., Neveu-Schwarz in CFT language) and periodic (Ramond) boundary conditions, respectively. This can be readily seen by comparing $T^{-1}$ and $\widetilde{T}$. These boundary conditions correspond to two replica-symmetric spin structures for the spacetime manifold. This is different from bosonic PT of fermionic systems [107, 110], where RN is given by sum over all possible spin structures. Essentially, the RNs associated with the two types of fermionic PT are identical to two terms in the expansion of bosonic PT in Ref. [107].

In what follows, we compute the two RNs for two partitioning schemes:

• •

Two adjacent intervals which is obtained by fusing the fields in $v_{A}$ and $u_{B}$. Hence, the RNs are given in terms of three-point correlators

[TABLE]

and

[TABLE]

where we introduce the fusion of unlike twist fields,

[TABLE]

• •

Bipartite geometry where the two intervals together form the entire system which is in the ground state. This time the RNs are obtained by further fusing the fields in $u_{A}$ and $v_{B}$ and the final expressions are therefore given by the two-point correlators

[TABLE]

and

[TABLE]

4 The spectrum of partial transpose

As mentioned, the first step to compute the tail distribution of the eigenvalues of partially transposed density matrix is to find its moments. To this end, it is more convenient to work in a new basis where the twist matrices are diagonal and decompose the partition function of multi-component field $\Psi$ to $n$ decoupled partition functions. For REE, this leads to $Z_{{\cal R}_{n}}=\prod_{k=-(n-1)/2}^{(n-1)/2}Z_{k,n}$, where

[TABLE]

The monodromy condition for the field around ${\cal T}_{k,n}$ and ${\cal T}_{k,n}^{-1}$ are given by $\psi_{k}\mapsto e^{\pm i2\pi k/n}\psi_{k}$. The calculation of the above partition function can be further simplified in terms of correlators of vertex operators using the bosonization technique in (1+1)d. For instance, in the case of REE, (86) can be evaluated by [114]

[TABLE]

where $V_{k}(x)=e^{-i\frac{k}{n}\phi(x)}$ is the vertex operator and the expectation values is understood on the ground state of the scalar-field theory ${\cal L}_{\phi}=\frac{1}{8\pi}\partial_{\mu}\phi\partial^{\mu}\phi$. The correlation function of the vertex operators is found by

[TABLE]

where $V_{e}(z)=e^{ie\phi(z)}$ is the vertex operator and $\sum_{j}e_{j}=0$. Hence, we can write for the partition function

[TABLE]

leading to the familiar result

[TABLE]

for the REE of 1d free fermions. Note that ellipses come from the proportionality constant in (89) which show sub-leading terms and may depend on microscopic details. In what follows, we apply the bosonization technique to evaluate ${Z}_{{\cal N}_{n}}^{({\rm ns})}$ and ${Z}_{{\cal N}_{n}}^{({\rm r})}$ similar to what we did for the REE. The scaling behavior of RNs in the lattice model is compared with the analytically predicted values of slopes (derived below) for various exponents $n=1,\cdots,7$ in Fig. 2, where the agreement is evident. We should note that the slope does not depend on the chemical potential $\mu$ in the Hamiltonian (9).

4.1 Spectrum of ${\rho^{T_{A}}}$

In the case of RN (69), we can carry out a similar momentum decomposition as

[TABLE]

where

[TABLE]

is the partition function in the presence of four twist fields. We then use (88) to compute the above correlator for various subsystem geometries. We should note that the following results only include the leading order term in the scaling limit, $\ell_{1},\ell_{2}\to\infty$, where $\ell_{1}$ and $\ell_{2}$ are the length of $A$ and $B$ subsystems, respectively.

4.1.1 Adjacent intervals

Here, we consider adjacent intervals (c.f. upper panel of Fig. 2(a)). The final result is given by

[TABLE]

where $f(x,y;q)=\frac{1}{2}\left[x^{2(q-1)(-2q+1)}y^{2q(-2q+1)}+x\leftrightarrow y\right]$. Notice that the exponents change discontinuously as a function of $k$. This can be understood as a consequence of the $2\pi$ ambiguity of the $U(1)$ phase that the Fermi field acquires as it goes around the twist fields. Essentially, we need to find the dominant term with the lowest scaling dimension in the mode expansion (see Appendix A for more details). Adding up the terms in the ${Z}_{{\cal N}_{n}}^{({\rm ns})}$ expansion, the final expression in the limit of two equal-length intervals $\ell_{1}=\ell_{2}$ is simplified into ${\cal N}_{n}^{({\rm ns})}=c_{n}\ln\ell+\cdots$ where

[TABLE]

where $N$ is a non-negative integer. It is worth recalling that for the bosonic systems, the spectrum of PT contains only positive and negative eigenvalues. As a result, we see even/odd effect for the moments. Here, however, the moments ${Z}_{n}^{(\text{ns})}$ have a cyclic behavior with a periodicity of six, which signals the possibility for the eigenvalues to appear with a multiple of $2\pi/6$ complex phase. As we will see below, this is indeed the case in our numerical calculations. We should also note that the above result can be obtained from the adjacent limit $v_{A}\to u_{B}$ of two disjoint intervals (92) as explained in Appendix B. Taking this limit is a bit tricky and was previously overlooked in Ref. [110], where it was incorrectly deduced that ${Z}_{{\cal N}_{n}}^{({\rm ns})}=0$ for two adjacent intervals.

We now discuss the spectrum of ${\rho^{T_{A}}}$ for two adjacent intervals. It is instructive to look at the many-body eigenvalues as obtained in (56) from the single-body eigenvalues of the covariance matrix (55). From the numerical observation that $\text{Im}(\nu_{k})\neq 0$, we may drop the $u_{l}$ factor in (56). Hence, the many-body spectrum simplifies to

[TABLE]

where

[TABLE]

and $\sigma_{k}=\pm$ is a sign factor. We should note that the complex and negative real eigenvalues come from product of $\omega_{I\sigma_{k}}$. This fact immediately implies that for every complex eigenvalue $\lambda_{j}$, $\lambda_{j}^{\ast}$ is also in the spectrum, since $\omega_{I-\sigma_{k}}=\omega_{I\sigma_{k}}^{\ast}$. Moreover, the negative eigenvalues are at least two-fold degenerate.

In the case of free fermions, we numerically observe that $\omega_{I\pm}\to|\omega_{I\pm}|e^{\pm i\frac{2\pi}{6}}$ as we go towards the thermodynamic limit $N_{A}=N_{B}\to\infty$. As a result, the many-body eigenvalues are divided into two groups: first, real positive eigenvalues, and second, the complex or negative eigenvalues which take a regular form $\lambda_{j}\approx|\lambda_{j}|e^{\pm i\frac{\pi}{3}s_{j}}$ where $s_{j}=1,2,3$. Figure 3(a) shows the numerical spectrum of ${\rho^{T_{A}}}$. To explicitly demonstrate the quantization of the complex phase of eigenvalues, we plot a histogram of the complex phase in Fig. 3(b) where sharp peaks at integer multiples of $\pi/3$ are evident. Due to this special structure of the eigenvalues, the moments of ${\rho^{T_{A}}}$ can be written as

[TABLE]

where $\{\lambda_{\alpha j}\},\alpha=0,1,2,3$ denote the eigenvalues along $\angle\lambda=\alpha\pi/3$ branches. Note that $\{\lambda_{0j}\},\{\lambda_{3j}\}$, i.e., positive and negative real eigenvalues, are treated separately, while $\{\lambda_{1j}\}$ and $\{\lambda_{2j}\}$ represent the eigenvalues for both $\angle\lambda=\pm\pi/3$ and $\angle\lambda=\pm 2\pi/3$ branches. A consequence of Eq. (4.1.1) is that there are four linearly independent combinations of the eigenvalues in ${Z}_{{\cal N}_{n}}^{({\rm ns})}$. This exactly matches the four possible scaling behaviors of ${Z}_{{\cal N}_{n}}^{({\rm ns})}$ from our continuum field theory calculations (100).

As a first characterization of the negativity spectrum, we compute the distribution of modulus of eigenvalues. To this end, it is sufficient to consider ${Z}_{{\cal N}_{n}}^{({\rm ns})}$ for multiples of $n=6N$ which is ${Z}_{{\cal N}_{n}}^{({\rm ns})}=\sum_{k}|\lambda_{k}|^{n}$. Substituting (100) for $b$ and $a$ in (27) and (30), we get

[TABLE]

where

[TABLE]

and $\lambda_{M}$ is the largest eigenvalue given by

[TABLE]

Figure 3(c) shows a good agreement between the analytical formula (104b) and the numerically obtained spectra for various subsystem sizes. We should note that there is no fitting parameter in (104b) and we only plug in $\lambda_{M}$ from numerics.

We can further derive the distribution of eigenvalues along different branches in Fig. 3(a). The idea is to analytically continue ${Z}_{{\cal N}_{n}}^{({\rm ns})}$ with $n=6N+m$ to arbitrary $n$ and solve the resulting four linearly independent equations generated by (4.1.1) to obtain the moments $\sum_{j}|\lambda_{s_{j}}|^{n}$ for each $s=0,\cdots,3$. This calculation relies on the assumption that $\lim_{n\to\infty}\frac{c_{n}}{n}$ does not depend on $m$, which is indeed the case in (100). Hence, we arrive at

[TABLE]

where $P_{\alpha}(\lambda)$ and $n_{\alpha}(\lambda)$, $\alpha=0,\cdots,3$ describe the distribution of eigenvalues along the $\angle\lambda=\alpha\pi/3$ branch. Here, $M$ and $\widetilde{M}$ encapsulate the coefficients

[TABLE]

$(a_{1},a_{2},\tilde{a}_{1},\tilde{a}_{2})=(\sqrt{\frac{3}{2}},1,\frac{1}{\sqrt{2}},\sqrt{3})$, and $\xi$ and $b$ are defined in Eqs. (105) and (106), respectively. Several comments regarding the phase-resolved distributions (107a) and (107b) are in order. The largest eigenvalue $\lambda_{M}>0$ is located on the real axis and hence only appears in $P_{0}(\lambda)$. The distribution of modulus is found by $(P_{0}+2P_{1}+2P_{2}+P_{3})$ which reproduces (104a). It is easy to check that the distribution is normalized and consistent with the identity $\text{Tr}{\rho^{T_{A}}}=1$,

[TABLE]

It is also possible to study the scaling of the maximum eigenvalue (in modulus) $|\lambda_{M}|$ along each branch. For the bosonic negativity, there are only two branches (positive and negative real axis) and it was found that the scaling of the maxima is the same in the thermodynamic limit [72]. In our case, for a given branch (labeled by $\alpha$) the maximum $|\lambda_{M}^{\alpha}|$ (with $|\lambda^{0}_{M}|\equiv\lambda_{M}$) can be extracted as

[TABLE]

where the result is independent of $\alpha$. This again implies the same scaling along each branch, up to a possible unknown constant due to non-universal coefficient that we are dropping in the above formulas (see Eq. (24)).

We compare the analytical results with the numerical simulations for each branch in Fig. 4. As expected, the numerical spectra reach the continuum field theory calculations as we make the system larger. We should point out that in contrast with the bosonic negativity spectrum and the entanglement spectrum which are given solely in terms of $I_{\alpha}(x)$, the modified Bessel function of the first kind, here the fermionic negativity spectrum contains the Bessel functions $J_{\alpha}(x)$ as well. Recall that unlike $I_{\alpha}(x)$ which is strictly positive for $x>0$, $J_{\alpha}(x)$ does oscillate between positive and negative values. Nevertheless, there is no issue in $P_{\alpha}(\lambda)$ which has to be non-negative, as the linear combinations of $I_{\alpha}$ and $J_{\alpha}$ in (107b) are such that they are strictly positive over their range of applicability within each branch.

4.1.2 Bipartite geometry

Here, we consider two intervals which make up the entire system as shown in the upper panel of Fig. 2(b). In this case, the branch points are identified pairwise as $u_{A}=v_{B}$ and $v_{A}=u_{B}$, where $\ell_{1}=v_{A}-u_{A}$. The partition functions in momentum space are found to be

[TABLE]

Similar to the adjacent intervals, the discontinuity in the $k$-dependence comes from the $2\pi$ ambiguity of the $U(1)$ monodromy (Appendix A). As a result, we have ${\cal N}_{n}^{({\rm ns})}=c_{n}\ln(\ell_{1})+\cdots$ where

[TABLE]

A benchmark of these expressions against the scaling of RN in numerical simulations is shown in Fig. 5(c). Because of the cyclic analyticity of the ${\cal N}_{n}^{({\rm ns})}$ modulo four, we expect to have the many-body eigenvalues along the real and imaginary axes. In other words, the complex phase of eigenvalues are multiples of $2\pi/4$.

We now derive the complex phase structure of many-body eigenvalues from the single particle spectrum. In the current case, the density matrix is pure leading to the identity $\gamma^{2}=\mathbb{I}$ for the covariance matrix. This property implies that the spectrum of the transformed covariance matrix (55) can be fully determined by the covariance matrix associated to the subsystem A, i.e., $\gamma_{AA}$ in Eq. (55). Hence, the single particle eigenvalues are given by

[TABLE]

and its Hermitian conjugate for $\nu_{k}^{\ast}$, where $\mu_{k}$’s ($k=0,\cdots,N_{A}$) denote the eigenvalues of $\gamma_{AA}$ [104]. Using (56), the many-body eigenvalues can be written as

[TABLE]

This decomposition has two types of factors: real positive and pure imaginary. Therefore, the many-body eigenvalues manifestly lie on the real and imaginary axes. Moreover, the many-body spectrum contains pairs of pure imaginary eigenvalues $\pm i\lambda_{j}$. The real negative eigenvalues are also two-fold degenerate since they are obtained from the product of even number of pure imaginary factors. In contrast, the real positive eigenvalues are not necessarily degenerate. As a result, the moments of ${\rho^{T_{A}}}$ take now the following form

[TABLE]

where $\{\lambda_{\alpha j}\},\alpha=0,1,2$ denote the eigenvalues along $\angle\lambda=\alpha\pi/2$. This expression in turn implies that there are three types of combinations of different branches for all $n$, which is again consistent with (121). By analytically continuing the three cases, we derive the moment $\sum_{j}|\lambda_{\alpha j}|^{n}$ for each branch. The resulting distributions are found to be

[TABLE]

where $M$ and $\widetilde{M}$ encode the coefficients

[TABLE]

$(a_{1},a_{2},\tilde{a})=(2,1,2\sqrt{2})$ and $\xi$ is defined in (105) with $b=-\ln\lambda_{M}=\frac{1}{6}\ln\ell$. As shown in Fig. 5, the above formulas are in decent agreement with numerical results.

Also in this case the maximum (in modulus) $|\lambda_{M}^{\alpha}|$ along the different branches can be evaluated through Eq. (114), giving (up to an unknown non-universal constant) $\ln|\lambda_{M}^{\alpha}|=-b$ independent of $\alpha$. Finally, also for the bipartite geometry, a consistency check is obtained from ${\rm Tr}\rho^{T_{A}}=1$, which simply follows from a calculation analogous to Eq. (4.1.1).

4.2 Spectrum of $\rho^{\widetilde{T}_{A}}$

Again, the first step to find the moments is the momentum decomposition of (71), yielding

[TABLE]

where the partition function

[TABLE]

is subject to modified monodromy conditions for the ${\widetilde{\cal T}}_{k,n}$ and ${\widetilde{\cal T}}_{k,n}^{-1}$, which are $\psi_{k}\mapsto e^{\pm i(2\pi k/n-\pi)}\psi_{k}$. This monodromy is different from the supersymmetric trace [118] (see Appendix C for the definition and more details).

4.2.1 Adjacent intervals

In this case, we find that

[TABLE]

It is important to note that for $k<0$, we modified the flux at $u_{1}$ and $v_{1}$ by inserting additional $2\pi$ and $-2\pi$ fluxes, respectively, where the scaling exponent takes its minimum value (c.f. Appendix A). Summing up $Z^{(\text{r})}_{k,n}$ terms, we get

[TABLE]

where

[TABLE]

for odd $n=n_{o}$, and

[TABLE]

for even $n=n_{e}$. As a consistency check, we show in Appendix B that the above formulae can be derived from two disjoints intervals as the distance between the intervals is taken to be zero. Notice that the even $n$ case is identical to the general CFT results [70]. Also, from (19) we arrive at the familiar result for the logarithmic negativity,

[TABLE]

For equal length intervals, we may write ${\cal N}_{n}^{({\rm r})}=c_{n}\ln\ell+\cdots$ where

[TABLE]

As expected for Hermitian operator $\rho^{\widetilde{T}_{A}}$, here the moments ${\cal N}_{n}^{({\rm r})}$ only depend on parity of $n$, i.e., whether $n$ is odd or even. This means that the eigenvalues are real positive or negative. We can also see this from the fact that the single particle spectrum is real. The many-body eigenvalues follow the form $\lambda_{\bm{\sigma}}=\prod_{\sigma_{k}=\pm}(1+\sigma_{k}\nu_{k})/2$, where $\nu_{k}$ are single-particle eigenvalues of the covariance matrix (59). As discussed in the previous section, we carry out the same procedure to derive the distribution from analytic continuation of moments (in this case there are only two branches). The final result reads

[TABLE]

where $\xi$ obeys the same form as Eq. (105) with a slight difference that $b=-\ln\lambda_{M}=\frac{1}{4}\ln\ell$. We present a comparison of the above expression with numerical spectrum of free fermions on the lattice of different lengths in Fig 6. There is a good agreement between analytical and numerical results.

We further find that, as it was the case for the bosonic negativity, the scaling of the minimum and maximum eigenvalue is the same. Finally, we confirm that the distribution probability is properly normalized such that $\int\lambda P(\lambda)d\lambda=\text{Tr}[\rho(-1)^{F_{A}}]$ and it is consistent with $\mathcal{E}=\frac{1}{4}\ln\ell$, Eq. (19), which follows from

[TABLE]

4.2.2 Bipartite geometry

In this case, we start by computing the correlator

[TABLE]

Here again, we have to minimize the scaling exponent for $k<0$ by inserting additional $2\pi$ fluxes (c.f. Appendix A). The RN is then found to be ${\cal N}_{n}^{({\rm r})}=c_{n}\ln(\ell_{1})+\cdots$ where

[TABLE]

From this, we derive the distribution of many-body eigenvalues to be

[TABLE]

where $\xi$ is given in (105) and $b=-\ln\lambda_{M}=\frac{1}{6}\ln\ell$.

We finish this part by a remark about the covariance matrix. Using the fact that $\gamma^{2}=\mathbb{I}$ for pure states, the covariance matrix (59) can be further simplified into

[TABLE]

Similar to the adjacent intervals, we can calculate the many-body spectrum out of eigenvalues of the above covariance matrix. We confirm that the numerical results and analytical expressions match. However, we avoid showing the plots here as they look quite similar to Fig. 6.

5 Conclusions

In summary, we study the distribution of the eigenvalues of partially transposed density matrices, aka the negativity spectrum, in free fermion chains. Taking the PT of fermionic density matrices is known to be a difficult task even for free fermions (or Gaussian states). However, recent studies [100, 101, 112] suggest that this difficulty could be circumvented if we use a different definition for partial transpose which is closely related to time-reversal transformation. In a matrix representation of a fermionic density matrix, e.g. in Fock space basis, the latter operation involves multiplying a $\mathbb{Z}_{4}$ complex phase factor in addition to the matrix transposition where the phase factor solely depends on the fermion-number parity of the state of subsystems to be exchanged in the transpose process. It turned out that the phase factor in the fermionic partial transpose lead to two types of partial transpose operation ${\rho^{T_{A}}}$ and $\rho^{\widetilde{T}_{A}}$. The difference is that ${\rho^{T_{A}}}$ is pseudo-Hermitian and may contain complex eigenvalues, while $\rho^{\widetilde{T}_{A}}$ is Hermitian and its eigenvalues are real. This is in contrast with the fact that the standard partial transpose ${\rho^{T_{A}}}$ is always a Hermitian operator which implies a real spectrum. In this paper, we presented analytical and numerical results for the negativity spectra using both types of fermionic partial transpose. In the case of $\rho^{\widetilde{T}_{A}}$, we find that the negativity spectra share a lot of similarities with those found in a previous CFT work [72]. However, in the case of ${\rho^{T_{A}}}$, we realize that the eigenvalues form a special pattern on the complex plane and fall on six branches with a quantized phase of $2\pi n/6$. The spectrum in the latter case is mirror symmetric with respect to the real axis, and there are four universal functions which describe the distributions along the six branches. The sixfold distribution of eigenvalues is not specific to complex fermion chain (described by the Dirac Hamiltonian) with $c=1$ and also appears in the critical Majorana chain with $c=1/2$. We further confirmed that our analytical expressions are applicable to the Majorana chain upon modifying the central charge $c$.

Given our free fermion results in one dimension, there are several avenues to pursue for future research. A natural extension is to explore possible structures in the negativity spectrum of free fermions in higher dimensions. It would also be interesting to understand the effect of disorder and spin-orbit coupling on this distribution. In particular, the random singlet phase (RSP) [119], which can be realized in the strongly disordered regime of one-dimensional free fermions, is characterized by logarithmic entanglement entropy [120, 121, 63] that is a hallmark of $(1+1)$d critical theories. An interesting question is how the negativity spectrum of critical RSP differs from the clean limit which was studied in this paper. Another direction could be studying strongly correlated fermion systems and specially interacting systems which have a description in terms of projected free fermions such as the Haldane-Shastry spin chain [122, 123]. Furthermore, it is worth investigating how thermal fluctuations affect the negativity spectrum in finite-temperature states. Finally, the negativity spectrum may be useful in studying the quench dynamics and shed light on thermalization.

Acknowledgments

The authors would like to acknowledge insightful discussions with Erik Tonni, David Huse, Zoltan Zimboras.

Funding information

SR and HS were supported in part by the National Science Foundation under Grant No. DMR-1455296, and under Grant No. NSF PHY-1748958. PC and PR acknowledge support from ERC under Consolidator grant number 771536 (NEMO). We all thank the Galileo Galilei Institute for Theoretical Physics for the hospitality and the INFN for partial support during the completion of this work. H.S. acknowledges the support from the ACRI fellowship (Italy) and the KITP graduate fellowship program. SR is supported by a Simons Investigator Grant from the Simons Foundation.

Appendix A Twist fields, bosonization, etc.

The Rényi entanglement entropy (REE) of a reduced density matrix $\rho$ is defined in Eq. (2). For non-interacting systems with conserved $U(1)$ charge, we can transform the trace formulas into a product of $n$ decoupled partition functions. Let us first illustrate this idea for REE [114]. We can diagonalize the twist matrix $T$ in Eq. (66) and rewrite the REE in terms of $n$-decoupled copies,

[TABLE]

where $\lambda_{k}=e^{i2\pi\frac{k}{n}}$ for $k=(n-1)/2,\cdots,(n-1)/2$ are eigenvalues of the twist matrix. In this new basis, the transformation rule $\Psi\to T\Psi$ for the field passing through the interval becomes a phase twist, i.e., $\psi_{k}\mapsto\lambda_{k}\psi_{k}$. Therefore, the REE can be decomposed into product of separate factors as

[TABLE]

where $Z_{k,n}$ is the partition function containing an interval with the twisting phase $2\pi k/n$. We reformulate the partition function in the presence of phase twisting intervals in terms of a theory subject to an external gauge field which is a pure gauge everywhere (except at the points $u_{i}$ and $v_{i}$ where it is vortex-like). This is obtained by a singular gauge transformation

[TABLE]

where $x_{0}$ is an arbitrary fixed point. Hence, for a subsystem made of $p$ intervals, $A=\bigcup_{i=1}^{p}[u_{i},v_{i}]$, we can absorb the boundary conditions across the intervals into an external gauge field and the resulting Lagrangian density becomes

[TABLE]

where the $U(1)$ flux is given by

[TABLE]

Note that there is an ambiguity in the flux strength, namely, $2\pi m$ (integer $m$) fluxes may be added to the right hand side of the above expression, while the monodromy for the fermion fields does not change. To preserve this symmetry (or redundancy), $Z_{k}$ must be written as a sum over all representations [124, 125, 126, 127]. The asymptotic behavior of each term in this expansion is a power law $\ell^{-\alpha_{m}}$ in thermodynamic limit (large (sub-)system size). Here, we are interested in the leading order term which corresponds to the smallest exponent $\alpha_{m}$.

As we will see in the case of entanglement negativity, we need to consider $m\neq 0$ for some values of $k$. Let us first discuss this expansion for a generic case. Let ${\cal S}_{n}$ be a partition function on a multi-sheet geometry (for either Rényi entropy or negativity). As mentioned, after diagonalizing the twist matrices ${\cal S}_{n}$ can be decomposed as

[TABLE]

where $Z_{k}$ is the partition function in the presence of $2p$ flux vortices at the two ends of $p$ intervals between $u_{2i-1}$ and $u_{2i}$, that is

[TABLE]

in which

[TABLE]

and $2\pi\nu_{k,i}$ is vorticity of gauge flux determined by the eigenvalues of the twist matrix. The total vorticity satisfies the neutrality condition $\sum_{i}\nu_{k,i}=0$ for a given $k$. In order to obtain the asymptotic behavior, one needs to take the sum over all the representations of $Z_{k}$ (i.e., flux vorticities mod $2\pi$),

[TABLE]

where $\{m_{i}\}$ is a set of integers and

[TABLE]

is the partition function for the following fluxes,

[TABLE]

and $\widetilde{\nu}_{k,i}=\nu_{k,i}+m_{i}$ are shifted flux vorticities. The neutrality condition requires $\sum_{i}m_{i}=0$. Using the bosonization technique, we obtain

[TABLE]

where $C_{\{m_{i}\}}$ is a constant depending on cutoff and microscopic details. We make use of the neutrality condition $-2\sum_{i<j}\widetilde{\nu}_{k,i}\widetilde{\nu}_{k,j}=\sum_{i}\widetilde{\nu}_{k,i}^{2}$ and rewrite

[TABLE]

where $\ell$ is a length scale. From this expansion, the leading order term in the limit $\ell\to\infty$ is clearly the one(s) which minimizes the quantity $\sum_{i}\widetilde{\nu}_{k,i}^{2}$. This is identical to the condition derived from the generalized Fisher-Hartwig conjecture [128, 129]. A careful determination of the leading order term for REE by a similar approach was previously discussed in Ref. [125, 126, 127].

We now carry out this process for ${Z}_{{\cal N}_{n}}^{({\rm ns})}$ in Eq. (92) for two adjacent intervals. Here, we need to minimize the quantity

[TABLE]

for a given $\nu=k/n=-(n-1)/2n,\cdots,(n-1)/2n$ by finding the integers $(m_{1},m_{2},m_{3})$ constrained by $\sum_{i}m_{i}=0$. For instance, let us compare $(0,0,0)$ with $(-1,1,0)$,

[TABLE]

So, we have

[TABLE]

Similarly, we find that

[TABLE]

In summary, we resolve the flux ambiguity by adding the triplet $(m_{1},m_{2},m_{3})$ as follows

[TABLE]

This leads us to write Eq. (95). Finally, similar derivation can be carried out to arrive at Eqs. (117), (132) and (144).

Appendix B Rényi negativity for disjoint intervals

In this appendix, we derive the RN associated with ${\rho^{T_{A}}}$ and $\rho^{\widetilde{T}_{A}}$ for two disjoint intervals and show that upon taking the distance between the intervals to zero, we recover the results for two adjacent intervals as discussed in the main text.

B.1 Moments of $\rho^{{T}_{A}}$

This geometry is characterized by $v_{A}-u_{A}=\ell_{1}$, $u_{B}-v_{A}=d$, and $v_{B}-u_{B}=\ell_{2}$ (c.f. Fig. 7(a)). The leading order term of the momentum decomposed partition function in the case of disjoint intervals is given by

[TABLE]

where

[TABLE]

Consequently, the RN is found to be

[TABLE]

We compare the above formula with the scaling behavior of the numerical results in Fig. 7(b), where we find that they match.

As a consistency check, we show that the RN between adjacent intervals can be derived as a limiting behavior of the disjoint intervals. However, we realize from (174) that $\lim_{d\to 0}Z_{k,n}^{(\text{ns})}=0$ (as is done also in Ref. [110]). A more careful treatment goes by considering higher order terms coming from different representations in (164)

[TABLE]

where $g(x,y;q)=x^{-2}\left(x/y\right)^{2|q|}+x\leftrightarrow y$ and $c_{ki}$ are coefficients dependent on the microscopic details. Next, we obtain the leading order term in the coincident limit $d=\varepsilon$, where $\varepsilon\ll\ell_{1},\ell_{2}$. To this end, we rewrite the above expansion (177) as

[TABLE]

where the scaling dimensions are

[TABLE]

As we see, for $|k/n|>1/3$, the second term is dominant. This immediately implies that upon taking $(\ell_{i}+d)\sim\ell_{i}$, we recover the original result (95).

B.2 Moments of $\rho^{\widetilde{T}_{A}}$

Similarly, we find the $k$-th contribution to the $n$-th moment of $\rho^{\widetilde{T}_{A}}$ to be

[TABLE]

which gives rise to the following form for the RN,

[TABLE]

where

[TABLE]

We compare the scaling behaviors of analytical expressions and numerical results in Fig. 7(c). As we see, they are in good agreement.

It is easy to verify that taking the adjacent limit $d=\varepsilon$ of two disjoint intervals in Eq. (181) leads to Eq. (133). We should note that in this case the leading order term in the momentum expansion (180) remains always the same (180) in contrast with the previous case (177).

Appendix C Partial transpose with supersymmetric trace

Let

[TABLE]

where $\widetilde{\text{Tr}}$ is the supersymmetric (susy) trace for the interval $A$. The susy trace is distinct from the regular trace in that the $T$ matrix which glues together $\rho^{\widetilde{T}_{A}}$ for fermions is given by (66), while the susy trace is similar to a bosonic trace (even though applied to a fermionic density matrix) where the $T$ matrix is given by (211) (see below). It is easy to see that $T^{n}=(-1)^{n-1}$ for the regular trace of fermionic density matrices whereas $T^{n}=1$ for the susy trace. Clearly, there is no difference between the susy and regular traces for even $n$ when considering $(\rho^{\widetilde{T}_{A}})^{n}$. The susy trace was used previously to define the susy entanglement entropies [118]. (See Refs. [130, 131, 132] for related works.) In terms of the partial transpose (13), the susy trace is simplified into

[TABLE]

which was studied by some of us in [112] and was shown to obey the same expressions as the bosonic negativity [70] for both even and odd values of $n$. In this appendix, we briefly report the results for various geometries. A technical point is that the monodromy of the field around ${\widetilde{\cal T}}_{k,n}$ for the susy trace is given by $\psi_{k}\mapsto e^{\pm i(2\pi k/n-\varphi_{n})}\psi_{k}$ where $\varphi_{n}=\pi$ or $\pi(n-1)/n$ for $n$ even or odd, respectively [110, 133].

1. Disjoint intervals: in this case the moments (189) become

[TABLE]

where $x$ is defined in (175), $c^{(1)}_{n}=-(n^{2}-1)/6n$, and

[TABLE]

2. Adjacent intervals: when the distance $d\to 0$, i.e., $x\to 0$ in the above expression, the moments take the form

[TABLE]

where

[TABLE]

for odd $n=n_{o}$, and

[TABLE]

for even $n=n_{e}$.

3. Bipartite geometry: finally, in this case one has

[TABLE]

where

[TABLE]

Appendix D Negativity of bosonic scalar field theory

As we have seen in the main text, calculating negativity boils down to computing correlators of twist fields. In this appendix, we briefly review the conformal weights of the twist fields in the complex scalar field theory,

[TABLE]

from which we can compute the correlators of twist fields and derive expressions for the entanglement of free bosons. Similar to fermions, the moments of density matrix in the coherent basis read as

[TABLE]

where

[TABLE]

For the moments of the partial transpose, we have

[TABLE]

In the case of free bosons, the moments can be written as a product of partition functions of decoupled modes,

[TABLE]

Hence, our objective here is to find the conformal weight of ${\cal T}_{k,n}$, ${\cal T}_{k,n}^{2}$, and their adjoints. It is worth noting that in the case of bosons $k$ takes positive integer values, $k=0,1,\cdots,n-1$. This is because the global boundary condition is periodic, i.e. the twist matrix obeys $T^{n}=1$. As usual in the conformal field theory, the computation goes by placing a twist field ${\cal T}_{k,n}$ at the origin which leads to a ground state ${\cal T}_{k,n}(0)\ket{0}$ where $\phi(z)$ and $\phi^{\ast}(z)$ are multivalued fields with the boundary conditions $\phi(e^{i2\pi}z)=e^{i2\pi k/n}\phi(z)$ and $\phi^{\ast}(e^{i2\pi}z)=e^{-i2\pi k/n}\phi^{\ast}(z)$. Next, we compute the correlator

[TABLE]

to find the expectation value of the energy-momentum tensor via

[TABLE]

Using the fact that

[TABLE]

we can read off the conformal weight $\Delta_{{\cal T}_{k,n}}$.

Let us start with ${\cal T}_{k,n}$ and ${\cal T}_{k,n}^{-1}$. The correlation function (215) can be directly computed by the mode expansion of $\phi$ field or can be simply derived from the asymptotic behavior $z\to w$ and $z\to 0$ or $w\to\infty$. The result is found to be [134, 135],

[TABLE]

which leads to

[TABLE]

We should note that doing this calculation for complex Dirac fermions, instead, leads to $\Delta_{{\cal T}_{k,n}}=k^{2}/2n^{2}$. So, the Rényi entropies are given by

[TABLE]

One can do a similar calculation for ${\cal T}_{k,n}^{2}$. In this case, the boundary condition is $\phi(e^{i2\pi}z)=e^{i4\pi k/n}\phi(z)$. For $k/n<1/2$, the result is identical to (219) up to replacing $k/n$ by $2k/n$. For $1/2<k/n<1$ however, the effective phase shift is $2\pi(2k/n-1)$ and we need to substitute $k/n$ in (219) by $2k/n-1$. This result can also be understood from the mode expansion. Consequently, we arrive at

[TABLE]

Using the following expression for the moments of partial transpose,

[TABLE]

we find

[TABLE]

and

[TABLE]

which are the familiar results [70].