Finite temperature negativity Hamiltonians of the massless Dirac fermion

The negativity Hamiltonian, defined as the logarithm of a partially transposed density matrix, provides an operatorial characterisation of mixed-state entanglement. However, so far, it has only been studied for the mixed-state density matrices corresponding to subsystems of globally pure states. Here, we consider as a genuine example of a mixed state the one-dimensional massless Dirac fermions in a system at finite temperature and size. As subsystems, we consider an arbitrary set of disjoint intervals. The structure of the corresponding negativity Hamiltonian resembles the one for the entanglement Hamiltonian in the same geometry: in addition to a local term proportional to the stress-energy tensor, each point is non-locally coupled to an infinite but discrete set of other points. However, when the lengths of the transposed and non-transposed intervals coincide, the structure remarkably simplifies and we retrieve the mild non-locality of the ground state negativity Hamiltonian. We also conjecture an exact expression for the negativity Hamiltonian associated to the twisted partial transpose, which is a Hermitian fermionic matrix. We finally obtain the continuum limit of both the local and bi-local operators from exact numerical computations in free-fermionic chains.


Introduction
In the last few decades, the study of entanglement turned out to be an optimal tool to investigate quantum field theories, quantum gravity models, condensed matter systems and quantum information theory [1][2][3][4][5][6]. Several entanglement measures have been studied in order to better probe the different features of a system. For example, in the context of pure states, the most important entanglement measures are the von Neumann and the Rényi entropies. They are defined as follows. Let us assume that our pure system is bipartite into A ∪ B and that the corresponding Hilbert space, H, factorises as H A ⊗ H B , where H A (H B ) is the Hilbert space containing the degrees of freedom of the subsystem A (B). The reduced density matrix of A is then obtained by tracing over the degrees of freedom of B and the entropies of ρ A , given by [7] S (n) ≡ 1 1 − n log Tr ρ n A , are good entanglement monotones on pure states. Despite the several successful applications of the entanglement entropy, they do not capture all of the properties of entanglement between the two subsystems. A more comprehensive measure of entanglement [8] is provided by the entanglement Hamiltonian K A , defined as the logarithm of the (appropriately normalised) reduced density matrix, i.e.
Being it an operator and not a scalar quantity, it is much more difficult to compute the entanglement Hamiltonian than the entropies, and it is known only for a limited number of cases. One of them is provided by the Bisognano-Wichmann theorem. This extremely general result applies to all Lorentz invariant quantum field theories in any dimension D + 1 and it states that, in the ground state, the entanglement Hamiltonian of the half-space A = x ∈ R D+1 x 1 > 0, x 0 = t = 0 is the generator of the Lorentz boosts preserving the Rindler wedge [9][10][11][12], i.e.
where T 00 is the energy density. This theorem provides a physical explanation of the Unruh effect [13] in terms of the entanglement of the vacuum. Although computing K A is a challenging task, we mention here that several achievements have been gained, for example in Gaussian states arising in lattice models of statistical mechanics and in quantum field theories (see e.g. [14][15][16][17][18][19]). If we shift our attention to mixed states, the entanglement entropies and Hamiltonian are not good entanglement measures, since they cannot distinguish between quantum and classical correlations. As a consequence, a plethora of alternative entanglement quantifiers have been proposed for dealing with the mixed states, even if in most cases these quantities are very difficult to compute even for few qubits (see e.g. [20,21]). A measure of entanglement in mixed states that has attracted a lot of interest is the logarithmic negativity [22][23][24], which is defined by doing a partial transpose operation on the reduced density matrix. To fix the ideas, let us consider a further bipartition of our subsystem A as A = A 1 ∪ A 2 and let the Hilbert space be factorised as H A = H A 1 ⊗ H A 2 . Then, if |e 1 i , |e 2 j are two arbitrary basis in, respectively, H A 1 and H A 2 , the partial transposition in, e.g., H A 1 acts as k , e 2 j |ρ A |e 1 i , e 2 l |e 1 i , e 2 j e 1 k , e 2 l | .
According to the Peres-Horodecki criterion [25,26], the presence of negative eigenvalues in the spectrum of ρ T 1 A is a sufficient, but not necessary, condition for the presence of entanglement between A 1 and A 2 . In light of this criterion, a quantifier of entanglement in mixed states is provided by the logarithmic negativity defined as [22][23][24] E ≡ log Tr ρ T 1 A .
Since the Peres-Horodecki criterion is only sufficient but not necessary, zero negativity does not imply the absence of entanglement. Nevertheless, this quantifier presents several advantages with respect to other entanglement monotones. On the one hand, from the quantum information point of view, as we have mentioned most measures of mixed state entanglement are very difficult to compute [20], while the negativity is computationally more efficient. On the other hand, from the field-theoretical perspective, the path integral construction of the partial transpose [27][28][29][30] makes it possible to calculate the negativity in field theories using the same conformal field theory (CFT) techniques adopted in the computation of the entanglement entropies.
In the framework of the mixed state entanglement, in [31] the negativity Hamiltonian has been introduced as the logarithm of the partially transposed reduced density matrix In particular, in [31,32] the negativity Hamiltonian of the ground state of the massless Dirac fermion has been explicitly computed for tripartite configurations, both in the presence and in the absence of boundaries, but a truly global mixed state has never been studied. In this work, we fill this gap computing the negativity Hamiltonian of several disjoint intervals at finite temperature, both on a finite size system and on the infinite line. This computation requires the knowledge of the entanglement Hamiltonian in the same geometry. Thus, Sec. 2 of the manuscript is devoted to a summary of the known results for the entanglement Hamiltonian of the massless Dirac field theory at finite temperature. This will set the stage for Sec. 3, where we will present the main results of this manuscript: we evaluate the negativity Hamiltonian of the Dirac fermion on the torus, i.e. at finite temperature on the circle, showing explicit results for a tripartite and a bipartite geometry. Our analytical findings will be checked in Sec. 4, where the continuum limit of the entanglement and negativity Hamiltonians is explicitly worked out starting from exact lattice numerical results. This allows us to check our predictions for the underlying Dirac field theory. In the same section we also propose an alternative definition for the negativity Hamiltonian corresponding to an Hermitian partial transpose reduced density matrix. After the discussion in section 5, we conclude the manuscript with two appendices containing computational details and additional results.

Finite temperature entanglement Hamiltonian
In order to derive the expression of the negativity Hamiltonian, we use the construction introduced in [31,32], where it was argued that the effect of the partial transposition amounts to exchange the extrema of the transposed interval in the expression of the entanglement Hamiltonian, taking into account that the fermionic field picks up an imaginary phase if it belongs to the transposed interval. Therefore, in this section, we present the known results for the finite temperature entanglement Hamiltonian of the free massless Dirac fermion in a multi-component region A, underlying the major differences with respect to the ground state.

Entanglement Hamiltonian on the torus
Let us consider a free massless Dirac fermion on a circle of circumference L at finite temperature 1/β, i.e., on a torus. In the imaginary time direction we impose anti-periodic (also called Neveu-Schwarz) boundary conditions, while in the spatial direction we choose either antiperiodic or periodic (Ramond) ones. Then, in a subsystem A = [a 1 , b 1 ] ∪ . . . ∪ [a n , b n ] composed of n intervals, the entanglement Hamiltonian is [33,34,37] where the signs + and − correspond, respectively, to the Ramond and Neveu-Schwarz sectors, p ∈ {0, . . . , n − 1}, k ∈ Z. The Hamiltonian (8) presents a local part, K loc A (β, L), proportional to the energy density T 00 , defined as (:: denotes normal ordering of the fields) with a weight given by the local entanglement temperature β loc (x) = 1/z (x), where [33,34,37] z(x; β, Here, = i b i − a i is the total lenght of the subsystem A and the additive constant term is only a shift which does not depend on x and, therefore, does not affect the expression for β loc (x), that we report explicitly In Eqs. (10) and (11), σ and ζ denote respectively Weierstrass' sigma and zeta functions and ϑ 1 is the Jacobi's elliptic theta function with nome q = e iπτ , τ = iβ/L (see Appendix B for their definitions). In particular, the expression in the first row of (10) is the result obtained in [33] while the one in the second row follows the conventions of Refs. [34,37]. While it is not obvious that the two alternative expressions coincide, one can show they are identical by using the properties of Weierstrass functions reported in Appendix B. In the rest of this manuscript we will adopt the conventions of [33] in terms of elliptic theta functions.
Regarding the non-local part K nl A (β, L) of Eq. (8), even in the case of one interval, this contains infinite terms proportional to the bi-local operator [18,35,36] In particular, the bi-local operator in (8) couples one point x with a single other conjugate pointx kp , given by the non-trivial solutions of the equations [33,34,37] z(x; β, L) − z(x kp ; β, L) + 2πk indexed by the integer k. One can see that for every fixed index k, Eq. (13) admits n solutions, indexed by p = 0, . . . , n − 1, each belonging to a different interval. In the following, we will use the index p = 0 to denote the solution of Eq. (13) such thatx k0 belongs to the same interval as x. With this convention, we see that for k = 0 Eq. (13) presents the trivial solution y =x 00 = x, which does not contribute to the non-local part K nl A (β, L) (see Eq. (8)). It is instructive to compare the entanglement Hamiltonian on the torus (8) with the one on the plane, i.e, of n intervals on the infinite line at zero temperature, given by [18,38] The local part of Eq. (14) is in form analogous to the one of Eq. (8), with entanglement temperature β loc (x) = 1/z (x) equal to the inverse of the derivative of the function The main qualitative difference of the Hamiltonian (8) on the toric space-time with respect to Eq. (14) is the structure of the non-local part K nl A . While in Eq. (8) the non-local part contains infinite terms, indexed by the integer k in Eq. (13), on the plane the bi-local part K bl A only contains n − 1 terms, calculated in the non-trivial solutions of the equation z(x) = z(x p ). In particular, for a single interval A = [0, ] the entanglement Hamiltonian (14) becomes completely local [18,19,[39][40][41] as expected since it is conformally equivalent to the Bisognano-Wichmann result (4). This shows that in general the entanglement Hamiltonian on the torus (8) is much more non-local than the analogous configuration on the plane [33,34,37].

Finite temperature entanglement Hamiltonian on the infinite line
We will now review the known results for the finite temperature entanglement Hamiltonian on the infinite line, i.e., on an infinite cylinder of circumference β in the time direction.
In [33,34,37], this Hamiltonian was obtained from the result on the torus (8) by taking the limit L → ∞.
Using the asymptotic expansion of the elliptic theta function ϑ 1 for q = e iπτ , τ = iβ/L → 0 (see Eq. (138) of the Appendix) in the expression (10) for the function z(x; β, L), we obtain where, using = i (b i − a i ), the contributions proportional to x cancel and we have introduced [33,34,37] z(x; β) = log The local term of Eq. (8) becomes proportional to the entanglement temperature [33,34,37] In the non-local component K nl A (β, L) of Eq. (8), instead, we can see that in this limit the denominator sinh(π(x −x kp + kL)/β) diverges for all k = 0 [33,34,37]. For this reason, the only conjugate points that contribute in this limit are the n − 1 non-trivial solutions of the equation [33,34,37] z( obtained as the limit of Eq. (13) with k = 0. This was expected by the fact that the cylinder can be conformally mapped into the plane, where the entanglement Hamiltonian is written in Eq. (14), which only contains n − 1 bi-local terms.
Putting all together, we find that the finite temperature entanglement Hamiltonian for a multi-component subsystem A = [a 1 , b 1 ] ∪ . . . ∪ [a n , b n ] on the infinite line is [33,34,37] with entanglement temperature β loc (x p ; β) given by Eq. (19). When specialising to a subsystem A made up of one interval, the entanglement Hamiltonian in Eq. (21) is purely local and in agreement with the result of [19,41], which reads As we mentioned earlier, since the cylinder is conformally equivalent to the plane, an alternative derivation of the finite temperature entanglement Hamiltonian on the infinite line in Eqs. (18), (21) consists in mapping the expressions (14), (15) on the plane to the cylinder. We find it worthwhile to also present this additional derivation as a non-trivial check of the correctness of Eq. (21) and because we will adapt a similar trick later in the manuscript. We first present how to map the entanglement Hamiltonian from the plane to a generic geometry and we later specialise this procedure to the cylinder. Let us consider a multi-component subsystem . . ∪ [a n , b n ] made up of n intervals in a geometry conformally isomorphic to the plane. In order to map it to the plane, it is convenient to switch to imaginary time w = x + it and consider, for simplicity, only the holomorphic component. Let then ξ(w) be the transformation from this geometry to the plane, with the subsystem A being mapped on the real line. On the complex plane, the holomorphic part of the entanglement Hamiltonian is given by the analytic continuation of Eq. (14) where the function z(w) = z(ξ(w)) is Eq. (15) evaluated in ξ(w), i.e.
with ξ(a i ), ξ(b i ) the extrema of the mapping ξ(A) of the subsystem A on the plane, and the conjugate pointsw p are the solutions of z(w) = z(w p ).
We first consider the mapping of the local part. Despite the fact that the holomorphic stressenergy tensor T is not a primary field, its transformation law only involves an additional function of w proportional to the Schwarzian derivative of ξ(w). When integrated in the entanglement Hamiltonian, this simply gives a constant factor which can be reabsorbed in the overall normalisation and can therefore be neglected. Considering also the Jacobian, the holomorphic part transforms as where we see that in the original geometry the entanglement temperature β loc (w) is given by the inverse of the derivative of Eq. (24) with respect to w. In order to find the transformation of the bi-local part, it is necessary to understand how the bi-local operator in Eq. (12) transforms under conformal mappings. In complex coordinates, the holomorphic bi-local operator takes the form Since the fermions ψ, ψ † are primary fields of conformal dimension ( 1 2 , 0), under the conformal mapping ξ(z) they transform as ψ(z) = ∂ξ ∂z 1/2 ψ(ξ(z)) (and analogously for the anti-holomorphic part). Replacing this transformation in Eq. (26) of the bi-local field, we find that in the original geometry it becomes Using the transformation of the bi-local operator in the entanglement Hamiltonian in Eq. (23), we obtain for the holomorphic part Putting together both the local and the bi-local components, we find that the holomorphic entanglement Hamiltonian in the original geometry takes the form and an analogous result can be also derived for the anti-holomorphic component. We remark that the expression of K A in Eq. (29) for multiple intervals under a generic conformal mapping to the complex plane is a novel result of this manuscript. In order to use the result of Eq. (29) for the finite temperature case, we recall that the cylinder is mapped into the plane under the transformation ξ(w) = e 2π β w , as shown in Fig. 1. In particular, at the time t = 0 in which we are interested in, the holomorphic and anti-holomorphic coordinate w coincides and the holomorphic and anti-holomorphic parts differ only in the operator. Substituting this mapping in Eq. (24), we reproduce the expression for z(x; β) at finite temperature reported in Eq. (18), which gives the entanglement temperature β loc (x; β) in Eq. (19). Regarding the bi-local part in Eq. (28), the weight function becomes which is also in agreement with Eq. (21), as expected. Therefore, we have used an alternative path to provide the results for the entanglement Hamiltonian of a disjoint set of intervals on the infinite cylinder. We stress that we find instructive to give this derivation here because we will use it also to evaluate the thermal twisted negativity Hamiltonian, which we introduce in the following section.

Finite temperature negativity Hamiltonian
In this section, we present the main analytical result of this paper, which is the field-theoretical prediction for the negativity Hamiltonian on a torus. After recalling the definition of the negativity Hamiltonian for fermionic systems, we review the construction introduced in [31,32] for its computation and we then extend it to the finite temperature case, showing two explicit examples. In particular, we find that in some cases the structure of the negativity Hamiltonian is more local than the one of the corresponding entanglement Hamiltonian.

General definitions
Let us consider a subsystem A = A 1 ∪ A 2 . As mentioned in the introduction, the negativity Hamiltonian N A in (7) is defined as the logarithm of the partially transposed reduced density matrix ρ T 1 A in (5), where we perform a transpose operation only in A 1 . The definition for such operation reported in Eq. (5) is appropriate for bosonic systems, but it turns out to be ill-suited for fermions: while the partial transposition of Gaussian bosonic states is still a Gaussian state, due to the anti-commutation relation, this is not the case for a fermionic one [42,43], and this makes the computation difficult even for Gaussian states [44][45][46]. For this reason, in [47][48][49][50][51] a more appropriate definition for fermionic systems has been introduced and the computational advantage one can gain is so remarkable that it has been employed in several contexts (see e.g. [52][53][54][55][56][57][58]). In order to motivate this definition, let us remark that, in a bosonic system, the partial transposition is equivalent to a partial time-reversal or a mirror reflection in phase space [26]. To see that this is indeed the case, we can consider a bosonic coherent state |α = e αa † |0 . On this state, the time-reversal transformation acts simply as the conjugation |α → |α * [26], therefore the relative density matrix goes into its own transpose For fermionic systems, the two transformations are not equivalent anymore. Under timereversal, a fermionic coherent state |ξ = e −ξc † |0 , ξ | = 0| e −cξ transforms as [47] |ξ ξ | −→ |iξ iξ| , which is different from the transposed density matrix because of the imaginary factor i. In light of this, one can define the partially time-reversed reduced density matrix ρ R 1 A , obtained by acting with Eq. (32) only in A 1 . This operation provides the fermionic logarithmic negativity E as although the spectrum of ρ R 1 A is not real in general [48]. An alternative definition for the fermionic partial transpose, called twisted partial transpose, has been introduced in [48] as where F A 1 = j∈A 1 n j is the number of fermions in the transposed subsystem A 1 . This new operator has only real eigenvalues and the logarithmic negativity given by [48] is a measure of the negativeness of the eigenvalues, exactly as for the bosonic partial transpose. In this sense, the twisted fermionic partial transpose has a more transparent interpretation of the fermionic negativity and allows for the measure of mixed-state entanglement also from its moments, in full analogy with the bosonic partial transpose [59,60].
Following the definition of negativity in Eq. (33), in [31] the fermionic negativity Hamiltonian N A has been defined as the logarithm of the partially time-reversed reduced density matrix ρ R 1 A , with an appropriate normalisation In order to compute this operator, in [31] it was introduced a physically motivated procedure to construct the negativity Hamiltonian (36) from the knowledge of the entanglement Hamiltonian, as we will review for the case of the ground state on the infinite line, i.e., on the plane. Moreover, in Appendix A we also apply the resolvent method to rigorously justify the construction of this operator. We can also define the twisted negativity Hamiltonian starting from Eq. (34) as but in this section we focus only on Eq. (36). We will come back to N A in Sec. 4. Before ending this section, we review the result for N A obtained in Ref. [31] (and we refer to the Appendix A for more details). Let us consider a multi-component subsystem A = [a 1 , b 1 ] ∪ . . . ∪ [a n , b n ] composed of n intervals and, to fix the ideas, let us reverse only one interval A 1 = [a j , b j ]; this case can be straightforwardly generalised to multiple reversed intervals. For this configuration, the entanglement Hamiltonian K A is given by Eqs. (14), (15) [18,38]. Under the path-integral construction of [27,28], the partial transpose has the net effect of applying a spatial reversal in the transposed interval. This can be understood in terms of CPT symmetry, since the time-reversal operation of Eqs. (31), (32) is equivalent to a parity transformation followed by a charge conjugation. This is implemented by exchanging the extrema a j , b j of the reversed interval in the expression of the entanglement Hamiltonian [31]. Under this procedure, the function (15) becomes Moreover, the effect of the partial transposition on the Dirac spinor ψ = ψ R ψ L is simply [32] (see also [47] and Appendix A). In the following section, we will show how this construction can be applied to compute the negativity Hamiltonian on the torus.

Negativity Hamiltonian on the torus
Starting from the result for the entanglement Hamiltonian in Eqs. (8), (10), by exchanging the extrema a j , b j we can obtain the negativity Hamiltonian on the torus. We remind that in the function z(x; β, L) in Eq. (10), it appears a term proportional to x and to the total length of the subsystem [33,34,37]. It is useful to write the subsystem length as = i (b i − a i ), since in order to obtain the correct negativity Hamiltonian it is necessary to exchange the endpoints of the reversed interval also in this expression. If we call 1 = j∈A 1 (b j − a j ) the total length of the partially reversed subsystem A 1 (for us, 1 = b j − a j ) and 2 = i∈A 2 (b i − a i ) the total length of A 2 , this procedure gives Analogously, Eq. (13), that determines the position of the conjugate points, becomes where again we have exchanged with 2 − 1 . For 1 = 2 , the negativity Hamiltonian on the torus has a non-local structure analogous to the one of the corresponding entanglement Hamiltonian in Eq. (8), containing infinite terms coupling different points where p ∈ {0, . . . , n − 1}, k ∈ Z and the function Θ 1 (x) is equal to 1 only for x ∈ A 1 , 0 otherwise In Eq. (41) we have introduced the "negativity temperature" and, analogously to Eq. (8), the signs + and − correspond respectively to the Ramond and to the Neveu-Schwarz sectors. On the other hand, when 1 = 2 , the dependence on the integer index k in Eq. (40) cancels out exactly and the solutionsx R kp with different k collapse on each another, giving a striking qualitative difference with respect to the entanglement Hamiltonian in Eq. (8). In Eq. (41), the bi-local terms with different k and same p are then calculated in the same conjugate pointx R p , leading to a bi-local structure with only n − 1 bi-local terms where we have introduced the (dimensionless) functions g R ± (z) defined by the infinite series In the Ramond sector (+ sign), Eq. (45) can be resummed to give where ψ q denotes the q-digamma function (see Appendix B for its definition), while in the Neveu-Schwarz sector (− sign) it reads To summarise, when the length of the reversed intervals is equal to the non-reversed one, the negativity Hamiltonian recovers a mild non-local structure given by a finite number of bi-local terms, while such a simplification does not arise in the entanglement Hamiltonian.
In the following we specialise to the case of n intervals lying on an infinite line at finite temperature (i.e. the space-time is a cylinder), and then we present explicit examples for the case of two intervals.

Finite temperature negativity Hamiltonian on the infinite line
The finite temperature negativity Hamiltonian on the infinite line can be obtained either by directly exchanging the extrema of the reversed interval in the related entanglement Hamiltonian reported in Eqs. (21), (18) or by taking the L → ∞ limit of the negativity Hamiltonian in Eq. (39), similarly to the limit reported in Eq. (17). By applying the exchanging procedure to Eq. (21), we find that the function z(x; β) in Eq. (18) reduces to and the n − 1 conjugate pointsx R p are found to be the non-trivial solutions of z R (x; β) = z R (x R p ; β). Thus, the finite temperature negativity Hamiltonian on the infinite line is where the negativity temperature β R loc (x; β) is given by and the bi-local terms are calculated in the n − 1 conjugate points obtained as the non-trivial solutions of z R (x; β) = z R (x R p ; β). As we also commented for the entanglement Hamiltonian, the negativity Hamiltonian only contains n − 1 bi-local terms.

Tripartite geometry
As a first explicit example regarding the negativity Hamiltonian on the torus, we consider a tripartite geometry made up of two intervals the length of A 1 and 2 = b 2 − a 2 the one of A 2 , and let us reverse the interval A 1 . Then, specialising Eq. (39) to this configuration we find while the conjugate point equation in Eq. (40) becomes We stress again that for 1 = 2 , the non-local structure of the negativity Hamiltonian drastically simplifies since the solutions of Eq. (52) do not depend on the index k, leading to a single bi-local term. We can now also consider some interesting limits of Eq. (51).
Finite temperature on the infinite line: If the two intervals which gives the negativity temperature The weight function of the bi-local operator reads As a further cross-check of our result, it is interesting to consider the zero-temperature limit β → ∞ of the negativity Hamiltonian. In this regime, we expect to retrieve the result for the tripartite configuration in the ground state, which was obtained in [31] by directly applying the exchanging procedure to the result in Eqs. (14), (15). Indeed, we see that taking the limit β → ∞ of z R (x; β) in Eq. (53), we reproduce the function on the plane found in [31] Regarding the conjugate points, again we see that they are given by the single non-trivial (57), finding the same conjugate point of [31] x as expected. For this geometry, i.e. two intervals on the plane, we provide a rigorous derivation of the result in Appendix A.

Bipartite geometry
We now study a bipartite geometry on the torus where A 1 = [0, 1 ] and the rest of the system is Notice that, differently from the case studied above, now the union A = A 1 ∪ A 2 of the reversed interval A 1 and A 2 is not a proper subset of the circle, but it covers all the system. Such a geometry can be obtained from the tripartite case of Sec. 3.4 by choosing a 1 = 0, b 1 = a 2 = 1 and b 2 = L. Taking this limits in the function z R in Eq. (51), we obtain where we have used the periodicity of the theta function ϑ 1 (z − π|q) = −ϑ 1 (z|q), while Eq. (52) for the conjugate points becomes The corresponding negativity temperature is provided by We again remark that for 1 = L/2 the dependence in k drops out from Eq. (60), and therefore all the infinite non-local solutions collapse into a single bi-local term with weight given by where g R ± are given in Eqs. (46) and (47), respectively. Let us stress that this represents an important result of this manuscript, since a bipartite system at finite temperature is a neat example of mixed state: in this case, the negativity is a genuine entanglement measure, differently from the entanglement entropy which mixes both quantum and thermal correlations. Therefore, the result for the negativity Hamiltonian provides the first operatorial characterisation of a thermal state. Let us now consider some interesting limits also for this bipartite geometry.
Finite temperature on the infinite line: Finding the theoretical prediction for the bipartite negativity Hamiltonian on the infinite line is more subtle than in the tripartite case of Sec. 3.4 because now A 1 and A 2 cover the full infinite line. The geometry of interest is and we reverse the interval A 1 . We can obtain this geometry from a three interval configuration on the infinite line , taking then the limit L → ∞ [51]. Specialising the function z R (x; β) in Eq. (48) to this geometry and taking the L → ∞ limit we find (up to x-independent terms) where now z R (x; β) reads This form differs from the one in Eq. (53) for the tripartite geometry, since now we find a term proportional to x. From this result we see that the negativity temperature is Since the geometry is made of three intervals, the equation for the conjugate points obtained from Eq. (48), with z R (x; β) = z R (y; β) is a polynomial of third order in y and one has the trivial solution y = x and also two non-trivial solutions y =x R ± , that in the limit L → ∞ read The bi-local inverse temperature corresponding to each conjugate pointx R ± is Another interesting limit we can study is when β → ∞, i.e. the zero temperature case, in which the state becomes pure. From Eq. (65), the negativity temperature is given by which is half of the weight function of the entanglement Hamiltonian for one single interval in the ground state in Eq. (16). The limit of Eq. (64) is and the two conjugate point in Eq. (66) arẽ In the limit β → ∞, the conjugate pointx R + (x R − ) diverges as O(β) for x < 1 /2 (x > 1 /2), and the bi-local operators calculated in this point do not contribute because the fermionic field ψ(x) vanish as x → ∞ [32,35]. In the other regions, instead,x R + andx R − are joined together to give the conjugate pointx R = x 1 /(2x − 1 ) in which the fermion does not vanish. Notice that, as expected, this conjugate point is precisely the only non-trivial solution of (69). We can explicitly compute the weight functions of the bi-local operators as As we can see, considering only the bi-local weights calculated in the region in which the conjugate points in Eq. (70) remain finite, the bipartite negativity Hamiltonian at zero temperature is We remark that, although one of the imaginary bi-local operators of the negativity Hamiltonian does not vanish, as β → ∞ the state becomes pure and [42,47]. As a consequence, we find The local part of the negativity Hamiltonian can be also rewritten as where I A 2 and I A 2 denote the identity operators on A 1 and A 2 , respectively, and are the entanglement Hamiltonians of the interval A 1 = [0, 1 ] (K A 1 ) and of its complement (K A 2 ). This result does not come as a surprise since a bipartite geometry at zero temperature is a pure state and one recovers that [47] Tr|ρ R 1 A | = Tr(ρ In other words, for a pure state the logarithmic negativity is equal to the Rényi entropy of order 1/2 defined in Eq. (2).

Numerical analysis
In this section we present exact numerical calculations on the lattice in order to compare them with our field-theoretical predictions. For Gaussian states, as those we are considering in this manuscript, it is possible to compute both the entanglement and the negativity Hamiltonian from the knowledge of the two-point correlation matrix restricted to the subsystem A, C A [15,17,[61][62][63][64][65][66]. For lattice fermions at finite temperature on the circle, C A is known both for periodic and for anti-periodic boundary conditions [67,68]. However, we stress that, as it has been studied in the literature [32,[69][70][71][72][73][74][75][76][77][78][79], comparing the lattice entanglement and negativity Hamiltonians with the analytical results is highly non-trivial. Indeed, while the terms of the field-theoretical predictions are localised around certain points (see, e.g. Eq. (14)), on the lattice the Hamiltonians are more delocalised, and this requires taking the continuum limit carefully. In this section, we first review how to recover the lattice entanglement and negativity Hamiltonians and we explain how to take the continuum limit of the lattice results. We then present numerical lattice computations, showing their good agreement with our predictions.

Lattice entanglement and negativity Hamiltonians for free fermions
On a circle of L sites, let us consider the tight-binding Hamiltonian where the lattice fermions satisfy the canonical anti-commutation relations and we impose either anti-periodic boundary conditions We can write down the Hamiltonian (77) in the Fourier modes c k , c † k and the dispersion relation of the tight-binding model (77) reads where the allowed momenta k depend on the boundary conditions, i.e., in the Neveu-Schwarz sector, the momenta are semi-integer while they are integer in the Ramond one Notice that, when L is divisible by 4, in the Ramond sector there are two zero-modes corresponding to the momenta k = ± L 4 . As discussed in [33,34,37,80,81], their presence is responsible for a non-local term in the ground state entanglement Hamiltonian. Choosing L = (2 mod 4) (i.e. divisible by 2 but not by 4), there are no zero-modes in the Ramond sector, while k = ±L/2 correspond to two zero-modes in the Neveu-Schwarz sector. To simplify the discussion, in the following we will focus on the case in which L is a multiple integer of 4. In terms of the energy ε(k) in Eq. (79), the two-point correlation matrix takes the form [51] Since the finite temperature state of free fermions is Gaussian, we can write its reduced density matrix in the subsystem A as [62] where h i,j plays the role of the matrix kernel of the entanglement Hamiltonian 2πK A . For Gaussian density matrices, h i,j is related to the two-point correlation matrix restricted to the subsystem A, C A , via Peschel's formula [15,17,[61][62][63][64] 1 where I A denotes the identity matrix in A and we have introduced the covariance matrix restricted to A By numerically computing Eq. (84) using the correlation matrices in Eq. (82), we obtain the lattice entanglement Hamiltonian in both the Neveu-Schwarz and Ramond sector. We stress that, due to the presence of the logarithm, the numerical computation of the formula (84) requires that the eigenvalues of C A are strictly included in (0, 1). For this reason, the numerical computation must be performed at high precision; in our study we used both the software Mathematica and the python library mpmath [82], keeping up to 300 digits. Peschel's formula (84) can be generalised to compute the negativity Hamiltonian, too. As we explained at the beginning of Sec. 3, the partial time-reversal ρ R 1 A of a Gaussian density matrix is still Gaussian, and this allows us to express the lattice negativity Hamiltonian as where now the kernel η i,j is non-hermitian. We consider a two-interval configuration , even though the generalisation to a multi-interval geometry is straightforward. Following [47], we write the blocks of the covariance matrix (85) as where Γ (σ,ζ) A i,j denotes i ∈ A σ , j ∈ A ζ . Under partial time-reversal in, e.g., the interval A 1 , the covariance matrix changes simply because of an imaginary factor for every index belonging to the reversed interval [47] The kernel of the negativity Hamiltonian η i,j in Eq. (86) is then related to the reversed covariance matrix Γ R 1 A in Eq. (88) in a way analogous to Peschel's formula (84) for the entanglement Hamiltonian [48] Since we are dealing with Gaussian states, also the relation between the twisted partial transpose ρ R 1 A in Eq. (34) and ρ R 1 A can be written as where the matrix U A = −I A 1 ⊕ I A 2 is related to the transformation (−1) F A 1 . Therefore, the kernel of the twisted negativity Hamiltonian, η, is given by As we mentioned, a proper comparison between the field-theoretical prediction and the lattice results requires a careful limiting procedure. To fix the ideas, we consider a configuration in which the field-theoretical entanglement Hamiltonian is local, such as a single interval on the plane in the ground state reported in Eq. (16) or at finite temperature in Eq. (22). The Hamiltonians in Eqs. (16) and (22) are completely local and proportional to the energy density T 00 . In light of this, one could naively expect that the corresponding entanglement Hamiltonian on the lattice would only have non-zero contribution on the first sub-diagonals h i,i+1 , h i+1,i . However, as studied for the first time in [69,70], this is not the case: On the lattice, higher hopping terms which couple fermions at longer distances are non-negligible and they must be included in order to recover the continuum limit. In the next subsection we review how to carry over this limit, showing why more care is needed for the case on the torus.

Continuum limit of the entanglement Hamiltonian
The key idea to take the continuum limit is to linearise the fluctuations of the lattice fermion c i around the Fermi points ±k F according to [70,83,84] where s is the lattice spacing (which will be put s = 1 in our numerical calculations) and we have introduced the continuous coordinate x = is. To lighten the notation, from now on, we redefine the momenta k introduced in Eq. (79) as k → k = 2k π/(Ls), such that k F = π/(2s) and we restore the lattice spacing s. The fields ψ L , ψ R are respectively the leftand right-moving components of a massless Dirac fermion, which describes the scaling limit of the tight-binding model in Eq. (77). Let us divide the entanglement Hamiltonian kernel h i,j in Eq. (83) in matrix blocks h (σ,ζ) i,j such that i ∈ A σ , j ∈ A ζ . In order to recover the local term of the entanglement Hamiltonian, one needs to consider the diagonal blocks h (σ,σ) . Following [70], one substitutes the linearisation in Eq. (92) in the expression of the lattice entanglement Hamiltonian in Eq. (83), obtaining where we expressed also the matrix element h (σ,σ) x,x+rs as a function of the continuous variable x. Since the massless Dirac field-theory presents conformal symmetry, one expects that in the continuum limit s → 0 the right-and left-moving fermions ψ R , ψ L will decouple. From Eq. (93) we can understand that the decoupling mechanism is due to the phases: the terms proportional to the product of left-and right-movers are multiplied by a strongly oscillating phase e ±ik F (2x+rs) and in the limit s → 0, these phases will average to zero, leading to the decoupling between ψ L and ψ R [70]. Dropping the highly oscillating terms and expanding in powers of s both the fields ψ L , ψ R and the matrix element h (σ,σ) x,x+rs we find We now plug the expansion (94) into Eq. (83) and we promote the sum over the index i to an integral over x. Integrating by parts the operator in the third row of (94), this term cancels out with the one proportional to the derivative of the matrix element ∂ x h in the second row.
In the second row, we recognise the number operator N (x) while in the last row the energy density T 00 (x) defined in Eq. (9). Thus, at leading order in the lattice spacing, we find that the diagonal blocks of the entanglement Hamiltonian can be written as [70] i where we have introduced the weighted sums over the matrix elements [70] S Let us compare Eq. (96) with the field-theoretical predictions for the entanglement Hamiltonian on the plane in Eq. (14) or at finite temperature in Eq. (21). Identifying the terms proportional to the energy density T 00 (x), in [70] it was verified that in the case of a single interval, the sum S loc (x) in Eq. (97) correctly reproduces the prediction for the local entanglement temperature β loc (x) both in the ground state (Eq. (16)) and at finite temperature (Eq. (22)). On the other hand, the term proportional to the number operator N (x) in Eq. (96) is expected to vanish at half-filling k F = π 2s . In this case, because of the particle-hole symmetry, the correlation matrix presents a checkerboard structure, inherited by the lattice entanglement Hamiltonian, which implies that Eq. (98) is identically zero. We are now interested in extending the analysis above to study free fermions on a torus, i.e. at finite temperature and size. The derivation of [70] reviewed in Eqs. (93), (94), (96) relies on the fact that all the matrix elements in the diagonal blocks contribute to the local term of the field-theoretical entanglement Hamiltonian (96). However, we have observed that the field-theoretical entanglement Hamiltonian in Eq. (8) contains infinite bi-local terms, even in the case of a single interval. This implies that summing over all matrix elements S loc (x) of Eq. (97) gives the wrong continuum limit, since we would be also including contributions that reproduce the bi-local terms of the entanglement Hamiltonian. It is therefore necessary to introduce a maximum cut-off R max in the sum in Eq. (97), to only include the local contributions. We show this in Fig. 2 for the local part of the entanglement Hamiltonian of one interval of length = 50 on the torus with L = 100 and β = 400. As we vary the cut-off R max ,  the agreement between the lattice bi-local weight in Eq. (97) and the theoretical prediction in Eq. (8) worsens. This non-local behaviour is also visible in Fig. 3, where we report the matrix plot of the entanglement Hamiltonian kernel h obtained via Eq. (84) for the case of a fermion on a torus at temperature β = 500 and system size L = 200 with anti-periodic boundary conditions. We see that besides the diagonal contributions, the matrix plot presents other terms located in the position of the conjugate points given by Eq. (13) for one interval.
In [75], the limit of the entanglement Hamiltonian was carried over also for the bi-local terms of a multi-interval entanglement Hamiltonian. Using Eq. (92) and again dropping the strongly oscillating contributions, we get In the second-to-last row we recognise the bi-local operator T bl (x, y) defined in Eq. (12), while the term in the last row is proportional to a different operator j bl (x, y) = j bl (x, y, 0) with which was already identified in [75]. In order to find the proper continuum limit, we now expand the field in position y around the conjugate pointx p , keeping only the term at leading order in s, obtaining [75] i j where we have again promoted the sum over the row index i to an integral over x and we have introduced the sums [75] If we now compare the limit of the off-diagonal blocks in Eq. (101) with the bi-local terms of the field-theoretical entanglement Hamiltonian (8), we see that the sum S bl (x) in Eq. (102) needs to reproduce the bi-local weight, since they are both proportional to the bi-local operator T bl (12). Analogously to what happens in the local case, we expect that the sum C bl (x) (103) multiplying the new operator j bl (x) (100) vanishes, since such an operator does not appear in the field-theoretical entanglement Hamiltonian (8). Also in the off-diagonal blocks, at half-filling k F = π 2s , the checkerboard structure of the lattice entanglement kernel h implies that Eq. (103) vanishes identically, simplifying the calculations.

Negativity Hamiltonian
In [32] it was argued that the limiting procedure of the lattice entanglement Hamiltonian h i,j reviewed in the previous section is almost identical to the one of the lattice negativity Hamiltonian η i,j of Eqs. (86), (89). Indeed, the limit only depends on the expansion of the This allows us to extract the negativity temperature from the lattice, in order to check our predictions of Sec. 3. The weight function of the local term can be read from Eq. (97), while the bi-local terms take different signs and imaginary factors in different intervals. In order to compare the continuum limit of the lattice negativity Hamiltonian, in the special case of two intervals, Eq. (102) must be modified as follows Also for the negativity, at half-filling Eqs. (98), (103) vanish identically. Now we can study the continuum limit of Eqs. (97) and (104) to check the field theory predictions for the negativity Hamiltonian, Eq. (41), for two disjoint intervals at finite temperature and size, in different regimes and both in a tripartite and bipartite geometry. In Fig. 4 we consider two adjacent intervals of equal length 1 = 2 , for several values of 1 and system size L and for different values of β, both with NS and R boundary conditions. In the left panel we find that the sum S loc over the higher hoppings is in perfect agreement with the field-theoretical local effective inverse temperature in Eq. (51). In the right panel, we report a similar analysis for the non-local term of the negativity Hamiltonian for the same geometry: we compare S bl in Eq. (104) with the field-theoretical weight function occurring in the bi-local term of the negativity Hamiltonian in Eq. (44), finding a good agreement. We stress that this geometry is quite interesting because the infinite non-local terms of the negativity Hamiltonian collapse on each other and we recover a bi-local structure, as we discussed in Sec. 3.2. This is also clear by studying the matrix plot of the kernel of the negativity Hamiltonian in   In Fig. 6, we consider again two intervals for different ratios of the length 2 / 1 = 0.5, 1, 1.5, with β/ 1 = 1/4. Here the system size is L = 20 1 , but since L β, this amounts to study a thermal tripartite geometry on the infinite line, whose analytical predictions are reported in Eq. (49). Indeed, both the left and the right panels confirm what we find analytically in Eqs. (54) and (56) for the local and bi-local terms of the negativity Hamiltonian, respectively. Before concluding the section, we want to check also the results for a bipartite geometry found in Sec. 3.5. In the top panels of Fig. 7, we consider a bipartition of a system of size L into two intervals of equal length, 1 = 2 = L/2, at inverse temperature β = L. This choice is particularly convenient because from Eq. (60) we can deduce that the infinite non-local terms are suppressed. Both the local and the bi-local component of the negativity Hamiltonian are in good agreement with Eq. (61) and Eq. (62), respectively. In the bottom panels, we consider a different geometry, and we perform a partial transpose operation with respect to the middle interval A 2 = [1, 1 ]. Since now A consists of three intervals, in the limit L → ∞, we have two conjugate pointsx R ± given by Eq. (66). We can find the continuum limit by studying We observe a good agreement with Eq. (65) for the local part (left) and Eq. (67) for the bi-local weight (right).

Twisted negativity Hamiltonian
While for the entanglement and negativity Hamiltonians we presented both known and novel field-theoretical predictions and we could compare them with the continuum limit of the lattice results, for the twisted negativity Hamiltonian defined in Eq. (37), there are no field theory results. To avoid confusion with the notation, we stress that we define the negativity Hamiltonian related to ρ R 1 A as N A and the one related to ρ R 1 A as N A . The advantage of studying ρ R 1 A is that it is an Hermitian operator, so the logarithmic negativity recovers its original meaning of measure of the negativeness of the eigenvalues. Although we do not manage to derive its form theoretically, we perform a numerical study on the lattice using the limiting procedure described in Sec. 4.2. This allows us to identify which operators appear in the continuum limit of the lattice twisted negativity Hamiltonian and we can formulate a conjecture for the local weight functions in the case of two identical intervals on the plane. We comment that this approach allows us to identify all the operators appearing in N A , contrarily to the analysis done in [31], where only the nearest neighbour negativity Hamiltonian has been considered.

Twisted negativity Hamiltonian on the plane
Let us first consider the twisted negativity Hamiltonian of the ground state on the infinite line, i.e, on the plane. The geometry under analysis consists of two adjacent intervals of identical length , and we perform a partial transpose operation on the first one, A 1 . As we did for N A , the continuum limit of N A is identical to the one of the entanglement Hamiltonian described in Sec. 4.2, since it depends only on the expansion of the lattice fermion in Eq. (92). However, differently from all the cases considered so far, we have numerically checked that even at half-filling k F = π 2s , the twisted negativity kernel η in Eq. (91) does not present a checkerboard structure. For this reason, also the terms proportional to the sums C loc (x) in Eq. (98) and C bl (x) in Eq. (103) have to be performed. This is the first difference with respect to Ref. [31], where the study of only the nearest neighbour terms prevented them from finding the operator C bl (x). This also confirms that, in order to recover the continuum limit correctly, a careful treatment of the long-range hoppings has to be taken into account. Therefore, besides the energy density T 00 (x) in Eq. (9) and the bi-local operator T bl (x, y) in Eq. (12), the continuum limit will contain also an imaginary local chemical potential term  proportional to the number operator N (x) in Eq. (95) and a term proportional to the operator j bl (x, y) defined in Eq. (100). Although we cannot derive the form of the weight functions of these operators explicitly, we provide a conjecture that very accurately matches numerical data on the lattice. Indeed, the twisted negativity Hamiltonian reads where the inverse negativity temperature β R loc (x) is given by the other weight functions are different and we report them herẽ The weight function of the number operator N (x) is the same that was conjectured in [31], while the weight functions for T bl (x,x R ) and j bl (x,x R ) are different and, we stress again, in order to recover them, it is important to sum over all the elements of the kernel of the negativity Hamiltonian, as done in Eq. (102). We also benchmark the analytical predictions from Eq. (106) in Fig. 8. The good agreement between the lattice computations and Eq. (106) supports our conjecture. Our prediction for equal intervals can be mapped into a geometry with adjacent intervals of different length using a Möbius transformation. for different values of the ratio 1 / 2 = 0.5, 1, 1.5. The system size is fixed as L/ 1 = 20 and we rescale the inverse temperature β such that β/ 1 = 1/4. The analytical predictions have been obtained by doing a conformal mapping from the plane to an infinite cylinder of circumference β in Eq. (113).

Möbius transformation
where we have used that the fermions ψ, ψ † transform as ψ(z) = ∂ξ ∂z 1/2 ψ(ξ(z)) (and analogously for the anti-holomorphic part). Therefore, taking into account Eq. (112) and the Jacobians of the transformation, we obtain the following expression for the twisted negativity Hamiltonian of two intervals of arbitrary length on the infinite line where β R loc (x) = 1/∂ x z R (ξ(x)) with z R given by Eq. (38). By doing another conformal mapping ξ(x) → e 2π β x in Eq. (113), we can obtain the result for two intervals on the infinite line at finite temperature, as shown in Fig. 1. We report a check of our conjecture in Fig. 9 for different ratios of the length 2 / 1 = 0.5, 1, 1.5, with β/ 1 = 1/4 and L = 20 1 . Beyond the good agreement, we observe that the weight function of the number operator N (x) drastically changes: the linear behaviour in x found at T = 0 becomes a kink interpolating from π for x < 1 to 0 for larger x. To summarise, starting from our conjecture for the twisted negativity Hamiltonian for two intervals of equal size on the infinite line, through a series of conformal mappings, we are able to find an expression also for the finite temperature case, which is a concrete example of a global mixed state.

Conclusions
In this manuscript we have continued the analysis initiated in Refs. [31,32] about the study of the negativity Hamiltonian, i.e. an operatorial characterisation of entanglement in mixed states. The most relevant novelty introduced here is the study of the entanglement in thermal states, which represent genuine examples of globally mixed states. Until now, the only configurations considered were non complementary subsystems at zero temperature. Here, we studied the negativity Hamiltonian of free massless Dirac fermions on a torus, for an arbitrary set of disjoint intervals at generic temperature. The structure of the negativity Hamiltonian exhibits a pattern similar to the entanglement Hamiltonian found in the same geometry in Ref. [33,34]: in addition to a local term, each point is non-locally coupled to an infinite but discrete set of other points. However, contrarily to what happens for the entanglement Hamiltonian, when the reversed and non-reversed subsystems have the same length, the bi-local solutions collapse on each other and we find only a finite number of bi-local terms, which couple each point only to another one in each other interval. We also analysed in detail the negativity Hamiltonian in a bipartite configuration. If the state is pure, the relation between the entanglement entropy and the negativity is well-known [27] and we retrieve it here. If the temperature is different from zero, a bipartite system is the first non-trivial example in which the negativity becomes essential to proper detect the quantum correlations. Also in this case, we found an infinite number of bi-local contributions, which reduce to one single bi-local solution only in the case of infinite system size. Our analytical findings are supported by exact numerical computations in a free-fermion chain. Another main result of this manuscript is the negativity Hamiltonian computed from the twisted partial transpose, cf. Eq (34). Through a careful numerical analysis, we identified the local and bi-local operators and their weight functions for two intervals on the infinite line both at zero and finite temperature. It would be interesting to derive analytically the conjectured formulae for the twisted negativity Hamiltonian, e.g. using the methods discussed in Appendix A. This study about the negativity Hamiltonian adds an important contribution to the operatorial characterisation of the mixed state entanglement, but there is still much work to do. For example, a challenging task is to exploit the mild non-locality of the negativity Hamiltonian together with the Hamiltonian reconstruction methods already used in [85][86][87] to reconstruct the negativity spectrum. Similarly, it is still an open problem to derive the conformal negativity spectrum [88] from the negativity Hamiltonian, as instead done for the entanglement spectrum in Ref. [89]. Another interesting direction is the study of the negativity Hamiltonian in higher dimensional systems, following what has been done for the entanglement Hamiltonian [74]. Finally, it would be also interesting to study whether one can define a notion of modular flow [11,90,91] for the partial transpose reduced density matrix and its eventual connections with the negativity Hamiltonian.

A The resolvent method for the negativity Hamiltonian
In [18], the field-theoretical prediction for the kernel H A of the entanglement Hamiltonian on the plane in Eq. (14) was obtained from the knowledge of the resolvent of the Green function C A restricted to the subsystem (see also [33,34,37,80,81,91]). In this appendix we show how to generalise the resolvent method of [18] to the negativity Hamiltonian in the case of multiple intervals on the plane, confirming the validity of the construction of Refs, [31,32] that we have used in Sec. 3.
For our purposes, we recast the resolvent method in terms of the partially reversed covariance matrix Γ R 1 A . To fix the ideas, we present the calculation for chiral fermions. Applying the partial reversal procedure in Eq. (88) to the Green function we find where the function Θ 1 (x), defined in Eq. (42), is equal to 1 only for x ∈ A 1 , 0 otherwise and P denotes Cauchy's principal value. Recall from the main text that the kernel of the negativity Hamiltonian can be related via Peschel's formula in Eq. (89) to the reversed covariance matrix Γ R 1 A . To apply Eq. (89) in the continuum theory, we first consider a single eigenvalue g of Γ R 1 A . For the entanglement Hamiltonian, in [18] it was used the fact that the spectrum of the Green function is real and contained in [0, 1]. In the case of the negativity Hamiltonian, we can use the knowledge that the eigenvalues of Γ R 1 A are contained in the unit complex disc |g| < 1 [48], as depicted in Fig. 10. Then, Peschel's formula for the single eigenvalue can be rewritten using Cauchy's theorem as where the branch cut of the logarithm is taken to go from −∞ to −1. Since |g| < 1, the contour of integration C in Eq. (115) can always be taken to avoid the branch cut (see Fig. 10) and therefore can be deformed continuously to integrate along the branch cut and on a small circle at infinity. Denoting the upper and lower branches of the complex logarithm as log + and log − respectively, and using the fact that the difference of the two branches is log + − log − = 2πi we find for every eigenvalue g of Since this holds for every eigenvalue, it holds also for the operator, leading finally to the expression for the kernel of the negativity Hamiltonian where we have introduced the resolvent of the partially reversed covariance matrix of Eq. (114) Note that throughout this appendix, 2πN A (x, y) corresponds to the continuum limit of η defined in Eq. (86). In order to find the explicit form of the resolvent in Eq. (118), we need to solve a singular integral equation. By construction, the resolvent (118) satisfies Multiplying both sides by (−)i Θ 1 (y) ζR(ζ; x, y) i Θ 1 (y) + (−1) Θ 1 (y) iπ P dz R(ζ; x, z) i Θ 1 (z) z − y = (−)i Θ 1 (y) δ(x − y) , we see that Eq. (120) has the form of a characteristic singular integral equation [92] a(y)φ(y) in the unknown function φ(y) = i Θ 1 (y) R(ζ; x, y), with the identification a = ζ, b(y) = (−1) Θ 1 (y) and f (y) = (−)i Θ 1 (y) δ(x − y). Comparing Eq. (120) with the analogous one for the entanglement Hamiltonian in [18], we see that the most important difference is the presence of the function b(y) = (−1) Θ 1 (y) in front of the Cauchy kernel, which changes sign if the interval is reversed. Now, we show that this function is precisely responsible for the inversion of the extrema a j , b j of the partially reversed intervals in the expression of the negativity Hamiltonian.
To solve Eq. (120), we introduce [92] G(y) = a(y) − b(y) a(y) + b(y) = ζ − (−1) Θ 1 (y) ζ + (−1) Θ 1 (y) = ζ − 1 ζ + 1 and the solution of Eq. (121) will be expressed in terms of the function [92] ω(y) = a 2 (y) − b 2 (y) exp 1 2πi P dz log G(z) z − y where z R is precisely the function in Eq. (38), obtained by exchanging the extrema a j , b j of the reversed intervals in the expression of Eq. (15). As we can see, the factor (−1) Θ 1 (z) in the second row of Eq. (123) is responsible for the exchange of the extrema in Eq. (38). The general solution of the characteristic singular equation (121) is [92] φ(y) = 1 a 2 (y) − b 2 (y) a(y)f (y) − b(y)ω(y) iπ P dz f (z) (z − y) ω(z) , which specialised to our Eq. (120) gives If we compare the resolvent for the negativity Hamiltonian on the plane in Eq. (120) with the one obtained in the context of the entanglement Hamiltonian in [18], we see that the main differences are the presence of the imaginary factors i Θ 1 (x) i Θ 1 (y) and the substitution of the function (15) with the one in Eq. (38) where the extrema of the reversed intervals are exchanged.
With the knowledge of the resolvent in Eq. (125), we can finally obtain the kernel of the negativity Hamiltonian by substituting it in Eq. (117). Changing variables as s = 1 2π log ζ−1 ζ+1 we find, formally (126) In the formal expression of the kernel N A (x, y), the Dirac delta is calculated in the solution of the equation z R (x) = z R (y). However, when dealing with the trivial solution y = x which corresponds to the local part of the kernel, Eq. (126) is proportional to the product of distributions δ(x − y)/(x − y) with coincident singular support. As discussed in [18], such an expression is ambiguous and it is necessary to regularise it. Following [18], the product is the distribution T that satisfies the algebraic distributional equation (x − y)T = δ(x − y), whose solution is T = −∂ x δ(x − y) + κ δ(x − y), where κ is an arbitrary constant which is fixed by requiring that the local part of N A is hermitian [18]. For this reason, we find it more convenient to explicitly antisymmetrise the kernel in the variables x and y, which cancels the κ δ(x − y) contribution. We also use the fact that the function z R in Eq. (38) has the property that it is monotonically decreasing in the reversed intervals A 1 and monotonically increasing outside, which implies for its derivative Then, by replacing Eq. (127) in the term of Eq. (126) corresponding to the trivial solution y = x we find which, when plugged in the expression for the negativity Hamiltonian reproduces the local part The n − 1 non-trivial solutions y =x R p of the equation z R (y) = z R (x) instead give rise to the bi-local terms. Explicitly anti-symmetrising the expression in the variables x, y gives leading to the bi-local part of the negativity Hamiltonian This resolvent procedure could be analogously extended to the case on the cylinder or on the torus considered in Sec. 3. Therefore, we can formally justify not only the construction introduced in [31] to compute the negativity Hamiltonian on the plane, but also at finite temperature or size, proving the correctness of the results found in this manuscript.

B Mathematical identities
We report here the main mathematical tools we have used throughout the manuscript. The Weierstrass zeta function is defined by [93] ζ It enters in the class of elliptic functions and it is quasiperiodic, i.e. it satisfies where P i , i = 1, 2, are the fundamental periods. In the case of interest for us, P 1 = L and P 2 = iβ. In order to prove the equality in Eq. (11), we have used the following representation of the Weierstrass zeta function through Jacobi functions For completeness, we also report here the definition of the Weierstrass sigma function used in Eq. (10) Also the equality in Eq. (10) can be proven by using the following property σ(x) = L π e ζ(L/2) x 2 L ϑ 1 π L x q ϑ 1 0 q .