Entanglement of a chiral fermion on the torus

In this paper we present the detailed calculation of a new modular Hamiltonian, namely that of a chiral fermion on a circle at non-zero temperature. We provide explicit results for an arbitrary collection of intervals, which we discuss at length by checking against known results in different limits and by computing the associated modular flow. We also compute the entanglement entropy, and we obtain a simple expression for it which appears to be more manageable than those already existing in the literature.


Introduction
In the last years, some ideas and results of quantum information theory have been applied in quantum field theory (QFT) with considerable success. For example, the strong subadditivity property of entanglement entropy has been used to show the irreversibility of the renormalization group flow in 2, 3 and 4 dimensions [1][2][3][4]. The interplay with quantum information is also prominent in holography, where the Ryu-Takayanagi formula [5] establishes a deep connection between entanglement and geometry.
Quantum information measures are based on the reduced density matrix, or equivalently on what is essentially its logarithm, the modular Hamiltonian. Among other applications, the knowledge of modular Hamiltonians was crucial for the proof of the averaged null energy condition [6], the derivation of quantum energy inequalities [7,8] and the formulation of a well-defined version of the Bekenstein bound [9]. Modular Hamiltonians also played a key role in applications to holography, the most notable case probably being the derivation of the linearized Einstein equations in the bulk from entanglement properties of the boundary CFT [10][11][12][13].
There are only few cases where modular Hamiltonians have been computed. For any QFT in the vacuum, the modular Hamiltonian of the half-space (or its causal development, the Rindler wedge) is proportional to the generator of boosts [14,15]. In the case of a CFT, one can use this result to derive the modular Hamiltonian of a ball by applying a special conformal transformation [16]. Further restricting to 1+1 dimensions, conformal maps from the plane to the cylinder can be used to obtain the modular Hamiltonian of an interval in other manifolds or states: a circle in the vacuum [17], the line at non-zero temperature [18]. Note that, since the boost generator is a local operator, all these modular Hamiltonians are local. In the case of free CFTs in the vacuum, some modular Hamiltonians are also known for multiple intervals: chiral fermion on the line [19] and the circle [20], chiral scalar on the line [21]. These modular Hamiltonians are not derivable from that of Rindler, and they turn out to be non-local. The non-locality disappears when restricting to a single interval, except in the case of the circle with Ramond boundary conditions, where the zero mode is responsible for the appearance of a non-local term even for a single interval.
In a recent paper [22] we computed a new modular Hamiltonian, namely that of a chiral fermion on a circle at non-zero temperature, i.e., on a torus in the Euclidean formalism (see [23] for a related result). The difficulty of the torus is that it cannot be mapped conformally to the plane. However, there is a very natural way to go from the torus to the plane: simply unfold the torus. This is, in broad terms, the idea of our method. More precisely, for free theories in Gaussian states there is a general expression for the modular Hamiltonian in terms of the resolvent of a certain two-point function, which is viewed as an operator [24]. Our strategy to compute the resolvent is based on the method of images applied to the calculation of the Euclidean propagator, which enables us to effectively unfold the torus and map the problem to a similar problem on the plane. In this paper we provide full detail on this calculation, and extend it by giving explicit results for multiple intervals. Interestingly, the modular Hamiltonian on the torus exhibits a completely non-local character, even for a single interval, regardless of the boundary conditions. We check the result by comparing it with known results in several limits, and we compute the modular flow, which inherits the non-locality of the modular Hamiltonian. Finally, we also use the resolvent to obtain a simple expression for the entanglement entropy of an arbitrary collection of intervals. The entanglement entropy on the torus had been previously computed in the form of a low-temperature and a high-temperature expansion [25,26]; we find a simple exact result which appears to be more suitable for applications.
The paper is organized as follows. In section 2 we present the expressions that give the modular Hamiltonian and the entanglement entropy in terms of the resolvent. In section 3 we compute the resolvent on the torus by the method of images. In section 4 we use the resolvent to compute the modular Hamiltonian, which we check in different limits against known results in the literature, and which we use to compute the modular flow. In section 5 we obtain a closed formula for the entanglement entropy and study its behavior, also checking the result in several limits. We have also included an appendix on the Weierstrass functions, which we use extensively throughout the paper.

Modular Hamiltonian and entropy from the resolvent
Consider a chiral fermion ψ on a circle of length L. The Hamiltonian is where the sign depends on the chirality. Suppose that the field is in a thermal state with inverse temperature β. The purpose of this paper is to compute the modular Hamiltonian H V corresponding to a subset V of the circle, which is related to the reduced density matrix ρ V by the equation In other words, H V is the Hamiltonian that the system should have in order for ρ V to be thermal with unit temperature. We will also compute the entanglement entropy Since the global state is Gaussian, correlation functions obey Wick's theorem. This is true everywhere in the circle, in particular within V , so the reduced density matrix is also Gaussian and hence the modular Hamiltonian has the form The result is [24,27] where both K V and G V are viewed as operators acting on functions on V . Note that G V is Hermitian and has its spectrum contained in the interval (0, 1), so K V is also Hermitian. Substituting (2.4) and (2.5) into (2.3) one obtains [24] The above two equations can be rewritten as This can be easily checked by explicitly performing the integrals in (2.7) and (2.8). Thus, the problem of computing the modular Hamiltonian and the entanglement entropy reduces to that of finding the resolvent of G V .

The method of images
Our strategy for computing the resolvent is based on the method of images applied to the calculation of the Euclidean propagator G, which is defined by for t − u ∈ (−β, β) and by analytic continuation for other values of t and u, where θ is the Heaviside step function. Depending on the spin structure chosen on the circle, the Euclidean propagator can be either periodic or antiperiodic in x with period L (these two alternative conditions are known respectively as Ramond and Neveu-Schwarz boundary conditions), and it is antiperiodic in t with period β. Due to these quasiperiodicity properties, we may view G as a section of a line bundle over a torus of circumferences L and β.
Taking the time derivative of the above equation and using (2.1) one obtains this equation says precisely that G(z, w) is analytic in z for z = w and has a simple pole at z = w with residue ±1/(2πi) (the latter part of this statement is shown by integrating (3.2) over a region containing (y, u) and then applying Green's theorem). In other words, where F is analytic in z for z − w ∈ (−L, L) × (−β, β). This equation is to be supplemented with the quasiperiodicity conditions for i = 1, 2, where P 1 = L, P 2 = iβ, ν 1 ∈ {0, 1} and ν 2 = 1. Eqs. (3.3) and (3.4) have a unique solution. Indeed, the difference ∆G between two solutions, viewed as a function of its first argument, is analytic for z − w ∈ (−L, L) × (−β, β) and satisfies (3.4), so it is analytic and bounded throughout the complex plane. By Liouville's theorem, such a function is necessarily a constant, so the antiperiodicity in the imaginary direction implies ∆G = 0. In order to find the solution, let us first look at the limiting case L, β → ∞, where the torus becomes a plane and the quasiperiodicity conditions (3.4) are replaced by the condition that G vanish at infinity. Since the only analytic function that vanishes at infinity is the zero function, the solution of Eq. (3.3) on the plane is  where Λ is the lattice and ν · λ = ν 1 λ 1 + ν 2 λ 2 = ν 1 λ 1 + λ 2 . Indeed, the function (3.6) clearly has the form (3.3), and one can easily check that it satisfies the quasiperiodicity conditions (3.4). This solution has to be taken as a formal one, because the above series is not absolutely convergent due to the slow decay of G 0 . The argument can be made rigorous by adding a mass m, which makes the Euclidean propagator on the plane an exponentially decaying function, and letting m → 0 at the end of the calculation.
Let now x, y ∈ V . It is clear from (3.1) that G V (x, y) = G(x, 0 + ; y, 0) (note that it is important to take the limit t → 0 from above, because if we take it from below we pick a delta function). With our identification of R 2 with C, we thus have so, by (3.6), where V Λ is a collection of segments distributed all over the complex plane, and The region V Λ is represented in Fig. 1 in the case where V is a single interval. The main reason why the method of images is useful for us is that an equation analogous to (3.9) holds for the powers of the operators involved, for any n ∈ N, where (A n )(u, v) denotes the kernel of the operator A n (not to be confused with [A(u, v)] n ). We can see this by induction. First, the above equation is satisfied for n = 1 (this is Eq. (3.9)). And second, if it holds for some n ∈ N we have which completes the proof. In the third equality we have used the translational invariance of G 0 , and in the fourth we have defined λ = λ + µ and z = z + µ. Eq. (3.12) implies that the method of images works for any function of G V which can be expressed as a power series. In particular, it works for the resolvent, (3.14) In the case of zero temperature, β → ∞, the terms with λ 2 = 0 do not contribute to the sum (3.6), so the lattice Λ effectively reduces to {mL, m ∈ Z} and, in consequence, the region V Λ reduces to an arrangement of segments in the real line. The resolvent R 0V Λ is well-known in that case [19], so we can use it to obtain R V via the above equation. In [22] we showed that, in fact, R 0V Λ can be easily computed for generic temperatures, where V Λ is a collection of segments distributed all over the complex plane. In the next section we review the argument in a slightly simplified form.

The resolvent for any set of segments in the plane
Let A be a collection of horizontal segments in the complex plane, see Fig. 2. The purpose of this section is to compute the resolvent of the operator G 0A with kernel for u, v ∈ A. As we will see, the computation is very simple: we will first obtain an expression for the powers of G 0A , then we will insert it into the power expansion of the x t Figure 2: A collection of horizontal segments in the complex plane.
resolvent and find that it is easy to perform the sum. For the square we have The second equality in (3.16) is immediately checked by performing the difference of fractions. Of course, the above integral is easily computed (it gives a logarithm), but for future purposes it is simplest to momentarily keep the integral expression in order to avoid dealing with the branches of the complex logarithm. Note that Indeed, if A +/− is a slight deformation of A that circumvents the point u ∈ A by taking the upper/lower path, we have where the contour in the last integral encircles u. Eq. (3.16), which gives the kernel of G 2 0A , can be rewritten as an operator equation, Now, using (3.18) and (3.20) it is a simple matter to check that the operator-valued function In turn, this implies for the n-th derivativeḠ (n) 0A = (−1) n n!Ḡ n+1 0A , as can be easily shown by induction. Noting thatḠ 0A (0) = G 0A , we thus obtain which is the desired expression for the powers of G 0A . Inserting it into the power expansion of the resolvent (which is a geometric series), we recognize the Taylor series ofḠ 0A , where in the last step we have used (3.18) to express ω ∓ A in terms of ω ± A and defined Finally, using (3.15) and the relation 1/( This is the resolvent for an arbitrary set A of horizontal segments in the complex plane. It coincides with the result of [19] in the case where A is contained in the real line.

The resolvent on the torus
Next we use the resolvent just computed, Eq. (3.25), to obtain the resolvent on the torus by the method of images, Eq. (3.14). From now on we will make extensive use of the Weierstrass functions, which are reviewed in appendix A (more extensive treatments are for example [28,29]). We suggest the unfamiliar reader to go through it before continuing.
where in the last step we have defined z = w − λ. The series on the last member of this equation is ambiguous because it is not absolutely convergent, as can be easily seen by the integral test. However, its second derivative is unambiguous, (3.27) where ℘ is the Weierstrass elliptic function. Therefore, where c and d are undetermined constants (which parameterize the ambiguity we found in (3.26)) and ζ is the Weierstrass zeta function. Note that property (3.18) and the fact that ζ has a simple pole with unit residue at the origin force c and d to be independent of the chirality. Let now N be the number of intervals in V , and let a i and b i be respectively the left and right endpoints of the i-th interval, that is, where σ is the Weierstrass sigma function (the iπ term comes from the interval containing x, which passes near a pole of the integrand). On the other hand, the quasiperiodicity of ζ, Eq. (A.9), implies where =´V dz is the total length of V and λ · ζ(P/2) = λ 1 ζ(P 1 /2) + λ 2 ζ(P 2 /2). Setting A = V Λ in the resolvent (3.25), substituting into (3.14) and using (3.28), (3.30) and (3.31), we obtain for the resolvent on the torus (3.25)), but the constant c remains. The latter is fixed by requiring the summand above not to diverge exponentially as |λ| grows. Indeed, taking into account that ζ is odd and satisfies ζ(z * ) = ζ * (z), so that ζ(L/2) is real and ζ(iβ/2) is imaginary, one easily sees that the only way to prevent the exponential divergence is to set so that the ambiguity in ω ± V Λ disappears from the resolvent on the torus. With this choice of c, the summand in (3.33) decays to zero as |λ| grows, but not quickly enough to make the series absolutely convergent, so there remains an ambiguity. However, note that, formally, F ± has the following properties: (i) it is analytic in z for z − w ∈ (−L, L) × (−β, β) except for a simple pole with residue ±1/(2πi) at z = w; and (ii) it is quasiperiodic, By the argument we gave under Eq. (3.4), there is only one function with these properties, so imposing them cures the ambiguity. The result is (recall that ν 2 = 1). Indeed, one can easily check from the properties of σ discussed in appendix A that the function (3.36) satisfies conditions (i) and (ii) above. To verify the latter condition one has to use the quasiperiodicity property which follows from the quasiperiodicity property of σ, Eq. (A.12), and the relation (A. 10) between the values of ζ at the half-periods. From (3.32) and (3.36) we finally obtain which is the resolvent on the torus and thus the main result of this section.

Modular Hamiltonian
We can now obtain the modular Hamiltonian by substituting the resolvent we just computed, Eq. (3.39), into (2.7). Note that the term with a delta function in (3.39) cancels from the modular Hamiltonian, because it changes sign under the transformation ξ → −ξ.
Changing the variable of integration from ξ to k(ξ) we find Surprisingly, the integral above is easily computed. To see this, note first from (3.38) that f is quasiperiodic in k, (in deriving the second equation in (4.2) we have also used the relation (A.10) between the values of ζ at the half-periods). Note that z is real because ζ(iβ/2) is imaginary. Note also that z(x) → ±∞ as x approaches the right/left endpoint of any interval of V , because σ vanishes at the origin. In appendix B we show that, in fact, z is a monotonically increasing function on each interval. Using the first of the above quasiperiodicity properties, we can rewrite the integral in (4.1) aŝ Except for a factor 2π, the sum in the second line is the Fourier series of the delta function on the circle, hence the last equality. Suppose now that x is fixed in the i-th interval, x ∈ (a i , b i ). Since z grows monotonically from −∞ to ∞ on each interval, for each n ∈ Z and j ∈ {1, . . . , N } the argument of the delta function has a unique zero y = x nj lying in the j-th interval. In other words, there is a unique x nj ∈ (a j , b j ) such that Clearly, x 0i = x. Noting from (4.1) that f (k; x, x) = 1, we may thus writê Let us compute the remaining integral. The quasiperiodicity properties (4.2), together with (4.5), imply and therefore, for the contour depicted in Fig. 3, we have  .7)). The integrand has only one pole inside the contour, which comes from a zero of the denominator in (4.1) and is located at k = −ν 1 β/(2 ) + iL/(2 ), as shown in the figure. Computing the contour integral by residues yieldŝ and inserting this result into (4.6) we obtain . (4.10) We can now substitute this expression into (4.1) to obtain the modular Hamiltonian. We find that the modular Hamiltonian is the sum of a local term and a non-local term, which come respectively from the first and second terms in (4.10). The local term is proportional to the product of distributions δ(y − x)/(x − y) and needs to be regularized. As explained in [19], the most general regularization gives δ , and the arbitrary function k can be fixed by requiring K V to be Hermitian. The result is . (4.14) Eqs. (4.11)-(4.13) give the modular Hamiltonian of a chiral fermion on the circle at non-zero temperature, for an arbitrary collection of intervals. Comparing the local term with the dynamical Hamiltonian (2.1), it is clear thatβ plays the role of a local inverse temperature. Note that the non-local term survives even in the case of a single interval, where the condition (n, j) = (0, i) reduces to n = 0. In that case, for x fixed, it has support in an infinite number of points, which tend to accumulate near the endpoints of the interval as shown in Fig. 4.

Limits
Let us now study some limits of the above result, to understand it better and to compare with known results in the literature. Consider first the limit L → ∞, with β finite. From Eqs. (A.13) and (A.14) we see that, in this limit, (4.3) becomes where we have dropped an irrelevant additive constant. On the other hand, in this limit the terms with n = 0 do not contribute to the sum in (4.13), so the non-local term takes the form where we have used the notation x j = x 0j . This gives the modular Hamiltonian of a chiral fermion on the line at non-zero temperature, for an arbitrary number of intervals. Note that the dependence on ν 1 , i.e., on the boundary conditions, drops out in this limit, as it should be. Further taking the limit β → ∞ of these expressions, one recovers the result of [19] for the modular Hamiltonian of chiral fermions on the line at zero temperature. The limit β → ∞, with L finite, is more subtle. From Eqs. (A.13) and (A.14) we see that, in this limit, (4.3) takes the form where again we have dropped an irrelevant additive constant. In order to properly take the limit of the non-local term, it is convenient to first apply it to a function f on V , Let us first compute the limit β → ∞ of the partial sum K (nl) V,1 f of terms with |n|L β. Noting from (4.5) that, for these terms, x nj = x 0j ≡ x j , we find where the infinite sums, which are not absolutely convergent, have to be taken in symmetric order because |n| is bounded by β. The first series above vanishes, and the second gives a trigonometric function that depends on ν 1 . The result is in agreement with [20]. What remains to be computed is the contribution K (nl) V,2 f to (4.18) from terms with |n|L ∼ β (the terms with |n|L β do not contribute in the limit β → ∞). Note first that, in this regime, the first two terms of the argument of the sinh in (4.18) can be dropped, because x − x nj is bounded by L and β is arbitrarily large. The remaining term can be rewritten using (4.5) as Note also, using (4.5) again, that z(x (n+1)j ) − z(x nj ) = 2π /β. In the limit β → ∞, this equation says that x (n+1)j is very close to x nj , so we can write Graphically, what happens for β large is that the curve plotted in Fig. 4 becomes arbitrarily steep, so the different solutions of (4.5) approach each other. Using these equations in (4.18) we obtain Now, since ∆x nj is vanishingly small in the limit β → ∞, in the Neveu-Schwarz case ν 1 = 1 this sum vanishes because consecutive terms cancel each other. Contrarily, in the Ramond case ν 1 = 0 the above sum becomes an integral, so we finally obtain , (4.25) in agreement with [20]. This emergence of continuous non-locality had already been noticed for the case of a single interval in [23]. Let us finally consider the case of a single interval, V = (a, b), and study the limit → L, i.e., the limit in which V becomes the whole circle, for arbitrary values of L and β. (where once again we have dropped an irrelevant additive constant), so the local temperature equals the true temperature,β = β. Note that, in this limit, z(x) does not diverge as x approaches the endpoints of the interval: the zeros of σ(a − x) and σ(b − x) cancel each other. In fact, the values of z at the endpoints are so small that Eq. (4.5) has no solutions inside the interval. What happens, as shown in Fig. 5, is that, as approaches L, the solutions of that equation approach the boundary of the interval, so that, in the limit, they have all reached the boundary. For x → a − , b + we have z (x) → ∞ and hencẽ β(x) → 0. It follows that the non-local term (4.13) vanishes in the limit → L, so the modular Hamiltonian becomes the dynamical Hamiltonian times β, as it should be.

Modular flow
Let us compute the modular evolution Ψ V of ψ, which is defined as i.e., analogously to the Heisenberg evolution with H V replacing H. For simplicity, we will restrict to the case where V is a single interval. The case of several intervals can be treated similarly. Differentiating the above equation with respect to s and using (4.11)-(4.13), we obtain where we have dropped the subscript j from x nj (the solution of Eq. (4.5) in the j-th interval) because there is only one interval. With this notation, recall that x 0 = x. The above is a non-local equation, but it can be rewritten as a system of coupled local equations by defining Ψ m (x, s) ≡ Ψ V (x m (x), s). Noting from (4.5) and (4.14) that which, after the change of variablesΨ m ≡ β (x m )Ψ m , simplifies to where M is the infinite-dimensional matrix with M mm = 0 and for m = n. Note that M is antisymmetric. It is to be regarded as a function of x, on which it depends via {x n }. Eq. (4.30) is a system of first-order local equations which can be solved by the method of characteristics. The solution is and P denotes path ordering. Note that, since M is antisymmetric, O ± is an orthogonal matrix. Note also that the curve x ± is well-defined because z grows monotonically from −∞ to ∞ and hence is invertible. Since the domain of z is the interval V , x ± (s) lies in V for all values of s. Undoing the change of variables and noting that Ψ V = Ψ 0 , we obtain .
as it should be. Let us now momentarily attach a superscript ± to each chirality to distinguish between them, and group them into the Dirac field ψ D = (ψ + , ψ − ). According to (4.34), the modular evolution Ψ DV of ψ D is .

(4.37)
Recalling that the Heisenberg evolution Ψ D of ψ D is given by Ψ D (x, t) = (ψ + (x+t), ψ − (x− t)), we can rewrite (4.36) as where y n (s) ± t n (s) = x ± n (s) and hence . (4.39) Note that (y n (s), t n (s)) is the intersection of two null lines which cross V at x + n and x − n , and hence it lies in D(V ), the causal development of V . Thus, according to (4.38), the modular evolution of ψ D (x) after a modular time s is a linear combination of the values of the Heisenberg Dirac field Ψ D at a countably infinite set of points in D(V ). As s changes, these points move, drawing an infinite set of curves which is represented in Figs. 6 and 7.

Entropy
The resolvent (3.39) can also be used to compute the entanglement entropy for the region V = N i=1 (a i , b i ) by substituting it into Eq. (2.8). In doing that, the terms proportional to the identity cancel, and changing the variable of integration from ξ to q = −k(ξ) one finds which is independent of the chirality. Since σ(z) ∼ z near the origin, the limit above is a derivative of the factor between square parentheses, where ζ ν (z) = (log σ ν ) (z) = ζ(z + L/2 + iν 1 β/2) − ζ(L/2) − ν 1 ζ(iβ/2). The left and right panels correspond respectively to Neveu-Schwarz (ν 1 = 1) and Ramond (ν 1 = 0) boundary conditions. At high temperatures the entropy is almost linear because thermal effects dominate over entanglement. At zero temperature, the Neveu-Schwarz entropy is symmetric with respect to the point = L/2 because the field is in a pure state. Contrarily, the Ramond entropy does not exhibit this behavior because the field is in a mixed state in which the zero mode has 1/2 probability of being empty or occupied. Correspondingly, the entropy in this case approaches the value S = log 2 as → L.
Note from the properties of ζ discussed in appendix A that ζ ν has no poles on the imaginary axis. Note also that it is odd, ζ ν (−z) = −ζ ν (z), and satisfies ζ ν (z * ) = ζ * ν (z). We can see that the ν-dependent term in (5.2) is independent of x, so the integral over V gives just the total length multiplied by the integrand ζ ν (iq ). For the other term, the integral in x can be performed introducing a cutoff to regulate the divergences coming from the entanglement across the endpoints of the intervals. The result is where we have defined .
Note that the dependence on the boundary conditions comes only through the second term in (5.4), which is real because ζ ν (iq ) is imaginary, and well-defined because ζ ν vanishes at the origin and diverges only linearly at infinity. Eqs. (5.4) and (5.5) give the entanglement entropy of a chiral fermion on a circle at non-zero temperature, for an arbitrary collection of intervals. This had been computed using the replica trick in [25,26] in the form of a low-temperature and a high-temperature expansion. We believe the above expression is much more practical for applications. In Fig. 8 we plot S V as a function of in the case of a single interval, for different temperatures which increase from bottom to top.
Let us now explore some limits of the above result. Consider first the limit L → ∞. From (A.13) and (A.14) we find that, in this limit, the entropy (5.4)-(5.5) takes the form which is independent of the boundary conditions, as it should be. Further taking the limit β → ∞ of the above expression, one recovers the result of [19] corresponding to the line and the vacuum state. The above expression can also be viewed as the first approximation to the entanglement entropy in the regime L β. Restricting it to the case of a single interval and letting → L yields as predicted by the Cardy formula. Consider now the limit β → ∞, with L finite. From (A.13) and (A.14) we find that, in this limit, (5.4)-(5.5) reduces to For Neveu-Schwarz boundary conditions (ν 1 = 1), the above formula implies S V = SV , whereV denotes the complement of V , as corresponds to a pure state which in this case is the vacuum. Contrarily, for Ramond boundary conditions (ν 1 = 0), the entropy does not have this symmetry property because g 0 (L − ) = g 0 ( ). This is not surprising: the chiral fermion with Ramond boundary conditions has a zero mode which does not contribute to the energy, and therefore its ground state is degenerate; it has degeneracy 2 because the zero mode can be either empty or occupied. Therefore, a Ramond chiral fermion at zero temperature is not in a pure state, but in a mixed state in which the zero mode has 1/2 probability of being empty or occupied. In fact, in the limit → L the above formula in the Ramond case gives S = 2πˆ∞ 0 dq tanh(πq) e 2πq − 1 = log 2, (5.10) in agreement with the above discussion. This is also shown in Fig. 8.
We also computed the mutual information between two diametrically opposed intervals of length 0 . This region and the corresponding mutual information as a function of 0 are shown in Fig. 9 (orange curve for Neveu-Schwartz and red curve for Ramond boundary  Figure 9: Behavior of the mutual information between the two diametrically opposed intervals shown in the figure as a function of their size 0 , for Neveu-Schwarz (orange) and Ramond (red) boundary conditions. In this plot L = 10 and β = 20.
conditions). The mutual information grows as the intervals get larger, which is consistent with the monotonicity property of mutual information. In the limit where the endpoints of the intervals coincide, the mutual information diverges as expected.

Conclusions
In this paper, we have presented the detailed calculation of a new modular Hamiltonian, namely that of a chiral fermion on a circle at non-zero temperature (i.e., on a torus in the Euclidean formalism), for an arbitrary collection of intervals. We have studied the properties of this modular Hamiltonian, checking that it reproduces known results in different limits and analyzing the corresponding modular flow. We have also obtained a simple exact expression for the entanglement entropy, which was previously known only in the form of a low-temperature and a high-temperature expansion [25,26]. Our approach is based on a novel and simple application of the method of images, which enables us to map the computation of the resolvent (the key object from which one can derive the modular Hamiltonian and the entanglement entropy) to a similar problem on the plane. As a result, we are able to obtain the resolvent by almost elementary calculations. We believe the method is very powerful and can be applied to other systems, especially free theories such as the scalar field, for which no known results at finite size and finite temperature exist in the literature.
Our result for the modular Hamiltonian is given in Eqs. (4.11)-(4.13). This is a highly non-local object, which couples each point of the region V to an infinite number of points within each interval of V , accumulating near the endpoints. The non-locality persists even in the case of a single interval, regardless of the boundary conditions, a behavior that had not been seen in any of the few modular Hamiltonians computed so far. The modular flow inherits this non-locality, and is described by an infinite set of curves which, as guaranteed by the Tomita-Takesaki theory, are contained in the causal development of V .
The entanglement entropy is given in Eqs. (5.4)-(5.5) for an arbitrary collection of intervals, and its behavior for a single interval is plotted in Fig. 8. In the Neveu-Schwarz case, the entropy of V becomes equal to that of its complement when the temperature tends to zero, consistently with the fact that, in that limit, the global state tends to a pure state (the vacuum). Interestingly, in the Ramond case the entropy does not exhibit this behavior, because the zero-temperature state is a mixed state in which the zero mode has 1/2 probability of being empty or occupied.
As mentioned above, we expect the applications of the method of images to go beyond the fermion case; we are currently following this strategy to study the entanglement properties of a scalar field on the torus. On the other hand, the resolvent we computed can find other applications beyond those presented here. For example, one can use it to study the effect of adding a small mass along the lines of [19]. Some features of the entanglement and Renyi entropies for a massive fermion on the torus have been studied before in [26], but we believe that more analytical progress can be achieved using the results of this paper.
with the lattice Λ = {λ 1 P 1 + λ 2 P 2 , λ i ∈ Z} is defined by where λ runs over Λ. One can easily check that ℘ is indeed elliptic with fundamental periods P 1 and P 2 . Within the fundamental period parallelogram, it has a double pole at the origin and no other poles, so it is of order 2. Note that it is even, ℘(−z) = ℘(z).
Together with periodicity, this implies that it is also symmetric, ℘(ω i − z) = ℘(ω i + z), with respect to the points ω 1 = P 1 /2, ω 2 = P 2 /2 and ω 3 = (P 1 + P 2 )/2, all lying in the fundamental period parallelogram. Since these three points are not poles, it follows that they are extrema of ℘. In fact, they are the only extrema of ℘ within the fundamental period parallelogram, because ℘ is an elliptic function of order 3 and hence has three zeros in each cell. Computing the Laurent expansion of ℘ around the origin and defining one finds that the function (℘ ) 2 − 4℘ 3 + g 2 ℘ + g 3 is analytic and vanishes at the origin. Since this function is elliptic, it must vanish everywhere, so that This is the Weierstrass differential equation. Note that it implies that ℘(ω 1 ), ℘(ω 2 ) and ℘(ω 3 ) are the three solutions of the cubic equation Q(t) = 0. In other words, for all t ∈ C. Equating the coefficients of each power of t, one obtains On the other hand, differentiating (A.4) with respect to z and using that ω 1 , ω 2 and ω 3 are necessarily first-order roots of ℘ , one finds that ℘(ω 1 ) ℘(ω 2 ) and ℘(ω 3 ) are first-order roots of Q, and hence they are all different from each other. Now, the case of interest for us is P 1 = L, P 2 = iβ. In this case the lattice Λ is invariant under complex conjugation, from which it follows that ℘(z * ) = ℘ * (z), (A.7) so ℘ is real when evaluated on the real line. Let us study the behavior of ℘ on the interval (0, L), which determines its behavior throughout the real line by periodicity. As we have seen, ℘ is symmetric with respect to the point ω 1 = L/2, where it has an extremum. This is the only extremum within the interval, because the other two extrema in the fundamental Figure 10: The Weierstrass ℘ function with P 1 = L, P 2 = iβ, evaluated on the real line.
period parallelogram, ω 2 and ω 3 , are not real. On the other hand, it is clear that ℘(x) → ∞ as x → 0, L, so x = L/2 is the absolute minimum of ℘ within the interval. Furthermore, ℘(L/2) > 0, so ℘ is positive. To see this, note first that, as a consequence of (A.7) and the periodicity of ℘, the three roots ℘(ω 1 ), ℘(ω 2 ) and ℘(ω 3 ) of the cubic polynomial Q are all real. Moreover, for x ∈ (0, L/2) we have ℘(x) ∈ (℘(ω 1 ), ∞) and ℘ (x) = 0, so, by (A.4), Q has no roots in the interval (℘(ω 1 ), ∞). Therefore, ℘(ω 1 ) is the largest root. Since the three roots are all different from each other, the first equation in (A.6) implies ℘(ω 1 ) > 0, as we wanted to show. The behavior of ℘ on the real line is shown in Fig. 10. The Weierstrass elliptic function has associated two quasiperiodic functions called the Weierstrass zeta and sigma functions. The Weierstrass zeta function ζ is defined by Clearly, this is an odd function satisfying ζ = −℘. The latter property, together with the periodicity of ℘, implies ζ(z + P i ) = ζ(z) + c for some constant c. Evaluating this equation at z = −P i /2 and using that ζ is odd yields c = 2ζ(P i /2), so that ζ(z + P i ) = ζ(z) + 2ζ(P i /2). (A.9) Note also that the poles of ζ are simple with unit residue, and they lie at the points congruent to the origin (i.e., which differ from the origin by an element of the lattice Λ).
Integrating ζ along the boundary of a cell and using the above quasiperiodicity property, one obtains P 2 ζ(P 1 /2) − P 1 ζ(P 2 /2) = iπ (A.10) provided that the ratio P 2 /P 1 has a pòsitive imaginary part (in the opposite case there is a minus sign on the right-hand side above). The Weierstrass sigma function σ is defined by (A.11) Figure 11: The Weierstrass ζ and σ functions with P 1 = L, P 2 = iβ, evaluated on the real line.
This is an odd function satisfying σ /σ = ζ. The latter property, together with (A.9), implies σ(z + P i ) = Ce 2ζ(P i /2)z σ(z) for some constant C. Evaluating this equation at z = −P i /2 and using that σ is odd yields C = −e ζ(P i /2)P i , so that σ(z + P i ) = −e ζ(P i /2)(2z+P i ) σ(z). (A.12) Note also that σ is analytic, has a simple zero at each point congruent to the origin and is non-vanishing everywhere else and satisfies σ (0) = 1. In the case of interest for us, P 1 = L, P 2 = iβ, the Weierstrass zeta and sigma functions satisfy f (z * ) = f * (z), like the Weierstrass elliptic function, so they are real when evaluated on the real line. We have ζ(x) → ∞ as x → 0 + and ζ(x) → −∞ as x → L − . Moreover, ζ is monotonically decreasing in the interval (0, L) due to the positivity of ℘. On the other hand, σ(x) = 0 at x = 0, L and σ(x) = 0 for x ∈ (0, L). Since σ (0) = 1, it follows that σ is positive within the interval (0, L). The behavior of ζ and σ on the real line is shown in Fig. 11. Let us finally study the behavior of the Weierstrass functions as one of the periods, say P j , goes to infinity while the other period, say P i , remains finite. In that case, all factors with λ j = 0 in (A.11) are equal to 1, so we can write σ(z) = z In the last step we have used the Euler product formula for the sine and recognized the Riemann zeta function ζ(2) = π 2 /6 in the exponent. From this expression we obtain ζ(z) = π P i cot πz P i + 1 3 (A.14) One can check that, in the case P 1 = L, P 2 = iβ, these limiting expressions have the qualitative features discussed above for generic values of L and β.

B Monotonicity of z
In the previous appendix we have shown that the Weierstrass zeta function is monotonically decreasing between any two consecutive poles in the real line. In fact, it also satisfies provided that x < y and there is no pole of ζ between x and y. This can be seen by using the relation between the Weierstrass elliptic function and the Jacobi theta function θ 1 , and also the expansion (log θ 1 ) (x) = cot x + 4 sin(2x) where q = e −πβ/L . Indeed, from the above expansion one sees that (log θ 1 ) (π/2) < 0, which, together with (B.2), implies ℘(L/2) > −2ζ(L/2)/L. Since x = L/2 is the absolute minimum of ℘, we have ℘(x) > −2ζ(L/2)/L for all x ∈ R, and integrating this relation yields (B.1). Let us now use this property to show that the function z defined in Eq. (4.3) is monotonically increasing on each interval of V , i.e., that z > 0. We have where we have used the relation (A.10) between the values of ζ at the half-periods. Rearranging the sum and using the quasiperiodicity of ζ, Eq. (A.9), we can rewrite f as

(B.5)
For each difference of zeta functions in the second line above, the arguments satisfy the conditions of (B.1). Using that property one immediately sees that f > 0, and hence z > 0, as we wanted to show.