Entanglement negativity in two-dimensional free lattice models

We study the scaling properties of the ground-state entanglement between finite subsystems of infinite two-dimensional free lattice models, as measured by the logarithmic negativity. For adjacent regions with a common boundary, we observe that the negativity follows a strict area law for a lattice of harmonic oscillators, whereas for fermionic hopping models the numerical results indicate a multiplicative logarithmic correction. In this latter case, we conjecture a formula for the prefactor of the area-law violating term, which is entirely determined by the geometries of the Fermi surface and the boundary between the subsystems. The conjecture is tested against numerical results and a good agreement is found.


I. INTRODUCTION
In recent years, ideas from quantum information theory have stimulated major developments in the field of strongly correlated systems. The entanglement properties of many-body states lies at the center of these studies. An important insight in this context is that, for ground states of local Hamiltonians, the entanglement between a subsystem and the rest of the system obeys an area law with a possible multiplicative logarithmic correction. 1,2 Moreover, the details of the groundstate entanglement scaling carries important information about the system, e.g., one can determine the universality class of one-dimensional critical models 3,4 or detect topological order. [5][6][7] While there has been numerous studies on the entanglement between a subsystem and its complement, much less is known about the entanglement between two regions that together do not constitute the entire system. The main reason is that the von Neumann entropy cannot be applied to study this case, since it is a measure of entanglement only for bipartite pure states. In the case of non-complementary regions, embedded in a larger system, one needs a different characterization, because the state reduced to the union of the subsystems is in general mixed. Among the various entanglement monotones for mixed states, entanglement negativity turns out to be a particularly useful measure. [8][9][10] It is easily computable for bosonic Gaussian states, 11,12 and recently also some results concerning fermionic Gaussian states have appeared. 13 Using these Gaussian methods, the tripartite ground-state entanglement negativity has been recently investigated for one-dimensional free bosonic [14][15][16] and fermionic models. 13,17,18 In this paper, we continue these surveys by considering two-dimensional lattices. When comparing the entanglement content of bosonic systems with that of fermionic ones, the dimensionality plays an important role. For critical one-dimensional systems, the entanglement between an interval and the rest of the chain scales logarith-mically with the length of the interval both for fermions and bosons. However, in higher dimensions, harmonic lattices (which can be viewed as free boson models) obey a strict area law, 19,20 while free fermions may violate the area law by a multiplicative logarithmic correction. [21][22][23][24] There is an appealing physical picture that gives an intuitive understanding of this difference. 25 For fermion models, the occupied and unoccupied modes in momentum space are separated by the Fermi-surface, characterized by a vanishing excitation gap. In the generic case, this is a (d−1)-dimensional surface of the d-dimensional Brillouin zone, around which the dispersion can be linearized. One can think that each patch of the Fermi surface is equivalent to a single gapless excitation, described by a 1+1 dimensional conformal field theory (CFT), which then leads to a multiplicative logarithmic correction to the area law. Indeed, the above argument implies exactly an entanglement entropy scaling that was already obtained based on the Widom conjecture. 22 In contrast, for harmonic lattices, the gapless bosonic modes are supported only on a single point (or on a discrete number of points) in the momentum space, which can give at most an additive logarithmic contribution to the entanglement entropy.
Here we set out to investigate the validity of the above intuitive physical picture also for the logarithmic negativity. We find that again a strict area law holds for harmonic oscillator systems, while in the case of free fermions our results indicate a multiplicative logarithmic correction. Moreover, using results obtained within a 1+1 dimensional CFT framework as an input, 14,15 we formulate a simple conjecture for the area-law violating term of the entanglement negativity, and present numerical evidence in its favor.
The paper is structured as follows. In Section II we recall the definition of logarithmic negativity and discuss how the partial transposition operation acts on bosonic and fermionic Gaussian states. CFT techniques concerning negativity are shortly reviewed in Section III, which are then used to calculate the negativity scaling for different subsystem geometries. The main results of the paper on two-dimensional models are presented in Section IV. We conclude in Section V with a short discussion of the results and their possible extensions. Various details of the analytical calculations are included in the two Appendices.

II. PARTIAL TRANSPOSE AND LOGARITHMIC NEGATIVITY
We will consider the ground state ρ of a many-body system defined on a lattice which is subdivided into three disjoint subsets A 1 , A 2 and B, such that A = A 1 ∪ A 2 and A ∪ B corresponds to the entire lattice. In case of a one-dimensional chain, two such tripartitionings are illustrated on Fig. 1. The reduced density matrix (RDM) of subsystem A is given by ρ A = Tr B ρ and its partial transpose w.r.t A 2 is defined through its matrix elements as In general, the result of the partial transposition is not a positive operator and the appearance of negative eigenvalues signals entanglement. 26,27 Based on this property, a suitable entanglement measure called logarithmic negativity 9 can be introduced as Despite being known as a computable measure of entanglement, the logarithmic negativity requires the knowledge of the full spectrum of the partial transpose ρ T2 A , which, in practice, is difficult to obtain for large many-body systems. A well-known exception is the class of Gaussian states. In this section we will review Gaussian techniques, which will be used later in Section IV.

A. Bosonic Gaussian states
Considering systems that are defined by a set of canonical coordinates {x n } and momenta {p n } with n = 1, . . . , N indexing the modes (or lattice sites), one can define continuous variable Gaussian states, also known as bosonic Gaussian states. Introducing the notation R 2n−1 = x n and R 2n = p n , bosonic Gaussian states are uniquely defined via their covariance matrix Γ kl = {R k , R l } , with higher order correlation functions factorizing according to Wick's theorem.
The reduced density matrix ρ A of a Gaussian state is again Gaussian and characterized by the reduced covariance matrix Γ A , where the indices are restricted to the subset A. Furthermore, the partial transposition has a particularly simple action on these states, since it can be represented as a partial time-reversal, flipping the sign of the momenta in the corresponding subsystem while leaving the coordinate variables unchanged. 11 In turn, the partial transpose ρ T2 A of the RDM is a Gaussian operator with covariance matrix 12 where is the diagonal matrix reversing the momenta in A 2 . Due to its Gaussianity, one has direct access to the full spectrum of the partial transpose via the symplectic spectrum {ν j } of Γ T2 A . In particular, the formula for the logarithmic negativity in Eq. (2) can directly be evaluated as This formula has been used in the earlier studies of negativity in various Gaussian many-body states. 28-31

B. Fermionic Gaussian states
Similarly to the bosonic case, the fermionic version of Gaussian states can also be defined, pertaining to a lattice system with creation c † n and annihilation operators c n satisfying canonical anticommutation relations {c † m , c n } = δ m,n . For a completely analogous treatment with the bosonic case, one can introduce Majorana operators a 2n−1 = c n + c † n and a 2n = i(c n − c † n ) and define the fermionic covariance matrix as Γ kl = [a k , a l ] /2. These two-point functions completely characterize a fermionic Gaussian state, as the higher-order correlations are given by the fermionic version of Wick's theorem.
Identically to its bosonic counterpart, the reduction of a fermionic Gaussian state to a subsystem A remains Gaussian with reduced covariance matrix Γ A . In sharp contrast, however, the partial transpose operation for fermions does not preserve Gaussianity. Nonetheless, it has been shown in Ref. [13], that in a suitable basis the partial transpose of a Gaussian RDM can be decomposed as the linear combination of two Gaussian operators. Indeed, the partial transposition with respect to A 2 leaves the modes in A 1 invariant and acts only on the ones in A 2 . Considering a product of n distinct Majorana operators M = a i1 a i2 · · · a in from subsystem A 2 , the transposition, in a particular basis, acts as M T2 = (−1) f (n) M where Using this definition, the partial transpose of a Gaussian RDM can be written in the form 13 Here O ± A are Gaussian operators with covariance matrices Γ ± A that are defined as where Although the spectra of O ± A can be constructed explicitly, the two operators do not commute in general and one has no direct access to the eigenvalues of ρ T2 A and, as a consequence, to the logarithmic negativity. Nevertheless, one can still extract some useful information from this form of the partial transposed RDM: the traces of its moments, i.e. Tr (ρ T2 A ) n . Indeed, factoring out Eq. (7), one is left with a sum of traces of products of Gaussian operators, each of which can be calculated explicitly. The steps of this procedure are summarized in Appendix A.
The moments of the partial transpose, despite not being entanglement measures, are the basic objects that are also attainable in CFT and were successfully used to characterize tripartite entanglement in 1D critical systems. [13][14][15] In this paper, we will also study quantities related to the moments to obtain an indication of the negativity scaling for 2D free fermion systems. To understand the role of the subsystem geometry in the 2D case, we first give a brief overview of the method employed in CFT to extract the moments of the partial transpose, and then apply it to simple 1D subsystem arrangements.

III. ENTANGLEMENT NEGATIVITY IN CFT
Conformal field theory provides a powerful machinery for the unraveling of universal properties of negativity scaling. Using CFT techniques, one could investigate the negativity in critical ground states, [14][15][16][17][18] in low temperature Gibbs states, 32,33 and even in non-equilibrium situations. 32,34-36 Moreover, some recent progress has been made in extending the technique to massive quantum field theories. 37 In this section we shortly review the main tools needed for such calculations, and use them to compare the negativity scaling for two different subsystem geometries depicted in Fig. 1. These results also constitute an essential input for our later studies of 2D free fermion models in Section IV.

A. The replica trick
The calculation of entanglement negativity in CFT relies on the path-integral representation of the partial transpose and on a clever application of the replica trick. 14,15 In the first step, one defines the ratio of the moments of the RDM and its partial transpose, as well as its logarithm Now, the crucial observation is that the ratio R n has a strong parity dependence due to the presence of negative eigenvalues of ρ T2 A . In particular, the trace norm in Eq.
(2) can be recovered by considering the series E ne on even integers n e and taking the limit We note that, instead of the ratio R n , one could simply use the moments of the partial transpose and their logarithms to obtain the entanglement negativity from the same limit as in (11). However, the definition of R n turns out to be more useful in various situations. On one hand, these ratios were shown to be universal in case of two non-adjacent intervals, i.e., R n depends only on the four-point ratio of the intervals through a universal scaling function corresponding to the given CFT. 15,38 On the other hand, while the same is not true for adjacent intervals, the ratios R n have a much clearer interpretation also in this case, as will be shown below.
In order to carry out the limit (11), one needs an explicit formula for E n and thus a method to calculate the traces in Eq. (10). This can be done by rewriting them as expectation values of products of twist fields T n and T n , permuting cyclically or anti-cyclically between the replicas. When both A i = [u i , v i ] with i = 1, 2 correspond to a single interval, one has 15 In other words, when considering moments of the RDM, the twist fields T n and T n have to be inserted at the start-and endpoints of the slits corresponding to A 1 and A 2 . For the partial transpose, the edges of the slit A 2 , and thus the corresponding twist field insertions have to be interchanged. Analogously, considering N non- , we can split them into two complementary sets I 1 and The moments are then given by where the partial swap operator S 2 acts as It should be mentioned that the general structure of these 2N -point functions of twist fields becomes rather involved for N > 2 and analytical results are only available for some special CFTs. 39

B. Adjacent intervals
We first consider the simplest situation with two adjacent intervals of lengths ℓ 1 and ℓ 2 , within an infinite onedimensional critical system. 15 One can then set u 1 = −ℓ 1 , v 2 = ℓ 2 and v 1 = u 2 = 0, hence R n can be written as and thus as the ratio of a three-point and a two-point function on the full complex plane. It is well known, that the twist fields T n and T n behave like primary operators with scaling dimension 4 In contrast, the numerator of Eq. (15) contains an insertion of a squared twist field T 2 n , whose scaling dimension shows a strong dependence on the number of replicas 15 Indeed, in case of n = n e even, the actions of T 2 n and T 2 n completely decouple the even and odd layers of replica sheets. Finally, one can use the CFT results for the two-point function where c n are non-universal constants. Similarly, the three-point function follows from conformal symmetry as where C TnT 2 n Tn are universal structure constants. Substituting into (15) and taking the logarithm, one arrives at which manifestly depends only on the scaling dimension ∆ (2) n . The entanglement negativity follows through the limit (11) as C. Embedded geometry In the case of an embedded geometry depicted in Fig. 1(b), one can perform a similar analysis as before. Choosing the subsystems as , the moments of the RDM and its partial transpose can be read off from Eq. (13) Due to global conformal invariance, the four point function appearing in (22) has the form 33,35 where F n is a scaling function depending only on the four-point ratio η defined as In writing Eq. (23) we have separated a term that becomes divergent in the limit η → 0, a behavior which can be found from the operator product expansion (OPE) technique. 33,35 Although the precise form of F n (η) is not known, such a definition ensures that it tends to constant values in both limits η → 0 and η → 1. We will consider a symmetric embedding, i.e., v 2 −u 2 = v 1 − u 1 = ℓ 1 and u 2 − v 1 = ℓ 2 , hence the four point-ratio is given by Dividing (23) by the two-point function, one obtains for the logarithmic ratio which, for η fixed and ℓ 1 , ℓ 2 ≫ 1, diverges logarithmically with a doubled prefactor compared to the result (20) for adjacent intervals. In particular for η → 0, i.e., for large separations ℓ 2 ≫ ℓ 1 between the two intervals of A 1 , one has F n (0) = C 2 TnT 2 n T n and thus the subleading term depends only on the structure constants of the corresponding OPE. Therefore, the entanglement negativity in the embedded geometry will be asymptotically twice as large as for adjacent intervals in Eq. (21), due to the two contact points between A 1 and A 2 . This is very reminiscent of the behavior of bipartite Rényi entropies with periodic vs. open boundary conditions. 4

IV. ENTANGLEMENT NEGATIVITY IN 2D
In the previous section it was recalled that for critical 1D systems the area law for negativity is violated logarithmically, and we also derived how the prefactor of this scaling depends on the subsystem geometry. In this section the analogous questions for 2D systems will be studied. In the bipartite case, the simple general connection between criticality and ground-state entanglement properties is lost. In fact, when considering the scaling of entanglement entropy of a subsystem, for many critical 2D systems, such as the harmonic lattice 23 or the Heisenberg model, 40 a strict area law holds. For other critical models, e.g., for free fermions 21,22,41 or interacting Fermi liquids, 42,43 multiplicative logarithmic corrections to the area law can still persist.
The above anomaly for fermionic models has its roots in the presence of a Fermi-surface, and its precursor can be traced back to 1D systems. Indeed, if the ground state of a free fermion chain is given by several disconnected Fermi seas instead of a single one, then the entanglement entropy between an interval of length ℓ and the rest of the chain gets multiplied by the number of Fermi points (boundaries of the Fermi seas). [44][45][46][47] Therefore, each Fermi point and subsystem boundary lends a 1/12 ln ℓ contribution of an independent chiral CFT to the entanglement. Hence the overall entropy is proportional both to the number of momentum-space boundaries of the Fermi seas and the real-space boundaries of the subsystem.
This argument can now be lifted to d-dimensional systems as follows. 25,48,49 Consider a ground state defined by the occupied fermionic modes whose border is given by a (d−1)-dimensional Fermi surface in momentum space, and suppose we are interested in the bipartite entanglement of a spatial region of linear extent ℓ. If each patch of the Fermi surface is considered as a source of chiral CFT excitations, with the direction of Fermi velocity given by the normal vector n q , then its entangling contribution along spatial direction n r should be proportional to |n q n r |. Summing up the contributions from the different patches of the Fermi surface ∂F and the real-space boundaries ∂A, one arrives at where the linear size of region A has been scaled out from the integral and the proper measure on the Fermi-surface has been taken into account.
Remarkably, the very simple argument leading to Eq. (27) gives the precise asymptotics of the area-law violating term in the entanglement entropy which has been tested numerically for a number of (2D or 3D) free-fermion systems 23,50-52 and recently even proved rigorously. 41 Moreover, it also accounts for the observation that the area law in 2D is restored whenever the Fermi surface degenerates to a number of points, 23,50 since ∂F is of zero measure in the integral (27) of the anomalous term. In fact, this is the very situation also for a harmonic lattice with short-ranged interactions, where the gapless bosonic modes are usually supported only on a finite number of points in momentum space. As a further proof of consistency, one should mention that some exotic Bose liquids, featuring a bosonic analogue of a Fermi surface, lead again to logarithmic violations of the area law. 53,54 Returning to the case of entanglement negativity, we expect that the same picture should be valid, which is supported by the numerics presented in this section. For harmonic lattices we find that a strict area law holds, while for fermionic systems the results are consistent with a multiplicative logarithmic correction, given by a formula that extends CFT results in the spirit of Eq. (27). To test its validity, we will both consider different subsystem geometries, as shown in Fig. 2, and different Fermi surfaces tuned by the anisotropy of the lattices.

A. Harmonic lattice
The Hamiltonian of a 2D lattice of coupled harmonic oscillators is given by where Ω 0 gives the strength of the harmonic confining potential at each site, whereas the oscillators are coupled through spring constants k m,n and the sum runs over all pairs m, n . Here we will restrict ourselves to the case of nearest neighbor couplings where the lattice sites are indexed by the pairs of integers R m = (i m , j m ). Introducing the notation q = (q x , q y ) for the wave vectors, the nonzero matrix elements Γ 2m−1,2n−1 = 2X m,n and Γ 2m,2n = 2P m,n of the covariance matrix read where the dispersion relation is given by The RDM of subsystem A has a reduced covariance matrix with nonzero elements given by 2X A and 2P A and the partial transposition acts only on P A → P T2 A by changing the signs of the momenta in A 2 . Due to the vanishing of cross-correlations between positions and momenta, the symplectic spectrum of the partial transposed covariance matrix is simply obtained by finding the eigenvalues {ν j } of matrix 2X A 2P T2 A . In turn, the logarithmic negativity can be evaluated through Eq. (5).
To be able to compare the results to the 1D case, especially to those obtained via CFT, we will be interested in a critical lattice system, i.e., in the limit Ω 0 → 0. Note, however, that one can not explicitly set Ω 0 = 0 due to a divergence in the matrix elements in Eq. (30) caused by the zero-mode of the lattice. We have thus used Ω 0 = 10 −3 in the calculations, and we observed that further decreasing Ω 0 has no visible effect on the results.
For simplicity, we consider only the vertical (Fig. 2a) and horizontal (Fig. 2b) partitions of a ℓ × ℓ square into two halves. The data for E are shown in Fig. 3 for the two different geometries, for different values of the vertical coupling k y and setting k x = 1. For the vertical partitioning (Fig. 3 left) and in the limit of uncoupled chains (k y = 0), one trivially recovers the c = 1 CFT result E/ℓ ∼ 1/4 ln ℓ + const. for the logarithmic negativity per area. However, already a small nonzero k y leads to a saturation of the curves, and thus to a strict area law of entanglement. This is indeed expected from the analogous result on the bipartite entanglement entropy, 23 which originates from the fact that there is only a single gapless mode within the Brillouin zone, i.e., one has Ω q = 0 only for q = (0, 0) for any k y = 0. Approaching the isotropic lattice k y → 1, the entanglement also becomes more evenly spread out in both directions, leading to a decrease (increase) of the logarithmic negativity across the vertical (horizontal) cut, as shown on the left (right) of Fig. 3.

B. Hopping model
The planar fermion hopping model is described by the Hamiltonian where, analogously to Eq. (29), we again consider nearest neighbor hopping only, with amplitudes t x and t y in the horizontal and vertical directions, respectively. Since the Hamiltonian given by Eq. (33) is particlenumber conserving, the problem simplifies considerably. Indeed, the basic quantities are the correlation functions C m,n = c † m c n which, in the ground state, are given by Here the integral goes over the Fermi sea, defined by q ∈ F if ω q < 0, with the single-particle dispersion obtained through diagonalizing Eq. (33) by a Fourier transform. Comparing to the formalism introduced in Sec. II B, one observes that the following relations hold for the elements of the covariance matrix Thus, in the presence of particle-number conservation, the spectrum of the reduced covariance matrix Γ A is directly related to that of the reduced matrix G A = 2C A − 1 1, which in turn determines the eigenvalues of the RDM ρ A . 55,56 In case of the partial transpose ρ T2 A , which can be written as a linear combination of two noncommuting Gaussian operators O ± A as in Eq. (7), it suffices to consider, instead of Γ ± A , the spectra of matrices G ± A to recover the eigenvalues of O ± A . Moreover, each moment Tr (ρ T2 A ) n can be obtained through determinant formulas involving only G ± A , as shown in Appendix A. Note that, since their linear size is half of the corresponding covariance matrices Γ ± A , their use is essential to reach large system sizes in our 2D calculations.
To further simplify the setting, we will consider only the geometries, depicted on Fig. 2, with A being a square of size ℓ × ℓ subdivided into rectangular subsystems A 1 and A 2 that share a common boundary. Before presenting our numerical results, and motivated by the results for the entanglement entropy for 2D free fermions in Eq. (27), we put forward a conjecture for the behavior of the logarithmic ratios where the geometric factor σ is given by Here the momentum and real space integrals have to be carried out along the Fermi surface ∂F and the common boundary ∂A 12 of subsystems A 1 and A 2 , respectively, and the linear size ℓ has been scaled out such that A becomes a unit square. Note that the prefactor (38) is chosen such that, in the limit t y = 0 of decoupled 1D chains, one recovers σ = c = 1 the central charge of the free fermion CFT, and Eq. (37) reproduces the result (20) for the adjacent intervals. Similarly, the doubling of the prefactor σ = 2 and hence Eq. (26) is recovered for the 1D embedded geometry. For the particular choices of the 2D partitions shown on Fig. 2 a) and b), a simple calculation (see Appendix B) yields whereas for cases c) and d) one trivially finds The numerical results for E n with n = 3 and n = 4, obtained by evaluating the determinant formulas in Appendix A, are shown in Figs. 4 and 5. In these figures we plotted −E n /ℓ on a logarithmic scale, to best visualize the expected behavior (37). The maximal linear size ℓ = 70 reached in our calculations is unfortunately still rather small to extract the slopes of the curves through fitting, as they become unstable due to the presence of subleading corrections. Instead, we simply plot the conjectured slope of the curves with dashed lines and compare it to the data sets.
In case of the vertical partitioning, one expects a slope which is independent of the vertical hopping amplitude t y . This is indeed nicely recovered from our data, shown in the left of Fig. 4, and even the numerical value of the slope fits very well to our conjecture. Decreasing t y , the area-law contribution increases and the curves are shifted upwards, which is the same trend observed for the negativity of the harmonic lattice, see left of Fig. 3. Additionally, the curves pick up oscillatory contributions with an increasing frequency.
On the other hand, for a horizontal partitioning the geometric prefactor σ b in (39) is a decreasing function of t y and thus the slope of the curves should tend to zero for t y → 0, i.e., when the partition becomes disconnected. This is clearly observed from our numerical data in the right of Fig. 4. In this case, however, it is somewhat more difficult to conclude about the correctness of our conjecture, as the data shows significant subleading corrections up to the reachable system sizes. Nevertheless, for n o = 3 one still finds a good agreement which, however, deteriorates for n e = 4. In fact, for 1D systems, the presence of unusual corrections whose magnitude increases with the index is well known for the Rényi entropies from CFT calculations 57 and was also observed for the moments of the partial transpose. 17 The data for the corner and the embedded partitionings are shown on the left and right of Fig. 5, respec-tively. Due to its geometric nature, the corner prefactor σ c is just the average of the vertical and horizontal ones, whereas the prefactor for the embedded geometry σ d is the double of the corner prefactor. This is indeed in very good agreement with the numerical data, especially for the embedded geometry which seems to be the least effected by subleading corrections.

V. CONCLUSIONS
We have studied the scaling of entanglement negativity in ground states of 2D free lattice systems between rectangular regions having a common boundary. While for harmonic oscillators a strict area law is obeyed, we found logarithmic corrections for the moments of the partial transpose in the fermionic case, which is completely analogous to the result for bipartite entanglement entropies. Based on this similarity and on CFT results for 1D systems, we conjectured a geometric form (38) for the prefactor governing the leading behavior (37) of the logarithmic ratios for the planar fermionic hopping model, and a comparison with numerical calculations shows a good agreement.
It would be interesting to find a strict proof for the form of the area-law violating term which, in the case of bipartite entanglement entropies, is related to the Widom conjecture 22 and has only been proved recently. 41 In contrast, for the moments of the partial transpose we do not even have an analogue of the method of Ref. [58] for 1D free fermions, where Rényi entropies are calculated us-ing the asymptotics of Toeplitz determinants. Although Tr (ρ T2 A ) n can also be cast as a sum of determinants, these are not of the Toeplitz type and we have not yet been able to find their asymptotics analytically. This would clearly be a necessary first step in order to understand the 2D results, and thus requires further studies.
There are a number of possible extensions of the setup presented in our work. Firstly, it would be interesting to see whether an interacting 2D Fermi liquid would show a similar negativity scaling as free fermions, which one would expect from the simple physical picture discussed in Section IV. The recent advances in numerical methods for evaluating entanglement negativity for interacting systems 59-63 cast some hope that this could be answered in the near future. Secondly, it would also be natural to investigate how the logarithmic correction to the negativity area law is rounded off when the Fermi surface degenerates to a number of points, and also to study corner effects. Thirdly, the question how the negativity decays with distance between non-adjacent subsystems should also be addressed. Finally, the possibility of detecting topological order via negativity 64,65 for freefermion systems, using the Gaussian toolbox presented here, is left for future study. Here we give the necessary formulae to evaluate the moments Tr (ρ T2 A ) n of the partial transpose in the ground state of a particle-number conserving free-fermion Hamiltonian. It will be additionally assumed that the correlation matrix C m,n = c † m c n is real, which holds for the models studied in the paper.
As described in the main text, the partial transpose can be given as a linear combination of two Gaussian operators, see Eq. (7). To simplify notation, we shall omit the subscripts A here and use O ± . In terms of the fermionic operators c † k and c k they are given by the quadratic form where σ = ± and Z σ ensures normalization, Tr O σ =1. The matrices in the exponent satisfy where G σ is defined through the correlation matrix C as and the subscripts refer to the reduction of matrices (rows and columns, respectively) to the corresponding subsystems A 1 and A 2 .
To obtain the n-th moment, one has first to factor out Eq. (7), which yields Note that the sum goes over all the possible assignments of {σ i } and i σ i is just the difference between the numbers of + and − terms in the corresponding factor. Using the fact that each O σi is a Gaussian operator given by Eq. (A1), one can apply determinant formulas for the traces of their products. Indeed, using the product relation for general Gaussian operators where exp(V ) = n i=1 exp(W σi ), and the trace formula 66 one obtains where the denominator is just the normalization factor n i=1 Z σi . Hence, the result (A4) is given by a weighted sum of 2 n determinants. In fact, it was recently pointed out that, for free fermions in 1D, each of these terms can be associated to partition functions on a higher genus Riemann surface with different spin structures. 18 Using the property O − = (O + ) † and thus the invariance of the formula under the exchange σ i → −σ i , the number of determinants to be evaluated can be reduced. Furthermore, using the relation (A2) between matrices W σ and G σ , each determinant in (A7) can be rewritten in terms of the latter ones. In particular, for n = 3 and n = 4 we obtain, after simple but lengthy algebra, the following expressions Tr (O 3 + ) = det Tr (O 4 + ) = det 1 1 + 6G 2 + + G 4 where Finally, one should note that very similar formulas can be applied in the general free-fermion case without particle-number conservation. 13 Then the matrices G σ should be exchanged with the modified covariance matrices Γ σ , given in Eq. (8), and one has to take the square root of the determinants. Note, however, that the resulting expressions then involve a sign ambiguity 67 which has its root in the underlying Pfaffian structure. 68

Appendix B: Calculation of the geometric prefactor
In this appendix we show how to calculate the geometric prefactor σ which determines the slope of the curves in Figs. 4 and 5, using a method very similar to the one presented in Ref. [52]. According to Eq. (38), the prefactor is given by a double integral over the Fermi surface ∂F and over the sur-face ∂A 12 separating subsystems A 1 and A 2 . The Fermi surface is defined through the single-particle dispersion (35) as q ∈ ∂F if ω q = 0 and is depicted on Fig. 6 for various anisotropies. Note that the Fermi surface is invariant under reflections q x → −q x and q y → −q y , it is thus enough to treat the first quadrant ∂F 1 in the integral (38) and multiply the result by four. Furthermore, for the geometries shown in Fig. 2 a) and b), the normal vector on the entire real-space surface ∂A 12 is constant and given by n r = (1, 0) and n r = (0, 1), respectively. Thus, in these cases the real-space integration can be dropped, whereas for the geometries in Fig. 2 c) and d) the results can be obtained trivially by combining those of a) and b).
The normal vector along the Fermi surface is given by and the path of integration can be parametrized as q x (θ), q y (θ), with the line element given by Furthermore, we can use the fact that the dispersion is constant (i.e. zero) along the Fermi surface which can be used to relate the line element in (B2) to the denominator of the normal vector in (B1). After proper cancellations, one finds the simple results Thus, the prefactors are simply related to the overall change in the components of the wavenumber as one sweeps along the Fermi surface. Looking at Fig. 6, the change in q y within the first quadrant is always given by π, independent of t y . On the other hand, the extent of the Fermi surface in the x-direction shrinks for larger anisotropies, with the locations of the endpoints given by q x = π/2 ± arcsin(t y ). Substituting into Eq. (B4), one obtains the results in (39).