The entanglement spectrum of chiral fermions on the torus

We determine the reduced density matrix of chiral fermions on the torus, for an arbitrary set of disjoint intervals and generic torus modulus. We find the resolvent, which yields the modular Hamiltonian in each spin sector. Together with a local term, it involves an infinite series of bi-local couplings, even for a single interval. These accumulate near the endpoints, where they become increasingly redshifted. Remarkably, in the presence of a zero mode, this set of points 'condenses' within the interval at low temperatures, yielding continuous non-locality.


INTRODUCTION
Amongst the predictions stemming from the interplay between Quantum Field Theory (QFT) and the causal structure of spacetime, one of the most robust is the celebrated Unruh effect: An accelerated observer in the vacuum measures a thermal bath, with a temperature proportional to its proper acceleration [1][2][3]. Intimately connected with the thermodynamics of black holes via Hawking radiation, this lies at the heart of our current understanding of the quantum nature of gravity [4]. Therefore, it is natural to explore its generalisations and investigate it further.
In recent years, these phenomena have been extended into the framework of quantum information theory. There, this temperature is understood as arising from the entanglement structure of the vacuum. Starting from a state ρ and some entangling subregion V , one defines the reduced density matrix ρ V by tracing out the complement of V . Then, just as the entanglement entropy S V = − Tr[ρ V log ρ V ] generalises the thermal entropy, the usual Hamiltonian is an instance of the more general concept of a modular (or entanglement) Hamiltonian K V defined via Originally introduced within algebraic QFT [5], the modular Hamiltonian has aroused much interest across a wide community due its close connection to quantum information measures. In the context of many body quantum systems, the spectrum of this operator is known as the "entanglement spectrum" and has been proposed as a fingerprint of topological order [6][7][8] and investigated in lattice models [9][10][11][12][13], as well as tensor networks [14][15][16]. In QFT, it is fundamental for the study of relative entropy [17,18] and its many applications to energy and information inequalities [19][20][21]. In the context of the AdS/CFT correspondence, it is instrumental in the program of reconstructing a gravitational bulk from the holographic data [22][23][24][25][26][27][28][29][30][31].
However, the modular Hamiltonian is known in only a handful of cases. The result is universal and local for the vacuum of any QFT reduced to Rindler space [3,32] and hence any CFT vacuum on the plane reduced to a ball [22]. For any CFT 2 , the same applies for a single interval, for the vacuum on the cylinder or a thermal state on the real line [33,34]. More generically, modular flows can be non-local, as is the case for multiple intervals in the vacuum of chiral fermions on the plane or the cylinder [35,36] and scalars on the plane [37]. The exact nature of the transit from locality to non-locality however is not fully understood, and remains an active topic of research.
In this paper we report progress regarding this problem, by providing a new entry to this list. We show that the chiral fermion on the torus (finite temperature on the circle) is a solvable model that undergoes such a transition between locality and non-locality. We compute the modular Hamiltonian by restating the problem as a singular integral equation, which in turn we solve via residue analysis.
Let us quickly quote our main result. For generic temperature, the modular Hamiltonian takes the form The local flow is of the standard Rindler form (7), with entanglement temperature given in (19). The novel result is the second term, given in (23) and depicted in Fig. 2, involving bi-local couplings between a discrete but infinite set of other points within the subregion. In the low temperature limit, the sector with a zero mode experiences a "condensation" of these points, resulting in a completely non-local flow.

THE RESOLVENT
We start by introducing the resolvent method, following [35,37,38]. For any spatial region V , the reduced arXiv:1905.05768v2 [hep-th] 15 Oct 2019 density matrix ρ V is defined as to reproduce expectation values of local observables supported on V . Now, for free fermions, Wick's theorem implies that it is sufficient that ρ V reproduces the equal-time Green's function for x, y ∈ V . This requirement fixes the modular Hamiltonian to be a quadratic operator given by [39] . This is specific for the free fermion. G V refers to the propagator as the kernel of an operator acting on functions with support on V .
As shown in [35] the modular Hamiltonian can be rewritten as in terms of the resolvent of the propagator, A derivation of (3) is provided in the supplemental material. In essence, it is the operator version of for a suitable choice of contour γ. In (4), the inverse of an operator is understood in the sense of a kernel, Thus, provided G of the global state and the entangling region V , this equation completely determines the resolvent R V and hence the modular Hamiltonian via (3).
To obtain the resolvent, let us first do the redefinition The convenience of this is that the first term of (5) will cancel the RHS of the previous equation, translating (4) into a singular integral equation All previous considerations hold for free fermions on a generic Riemann surface. The simplest case is the plane where the solution of (6) is a standard result [40], which was used by [35] to derive the corresponding modular Hamiltonian. They found that for multiple intervals, it consists of a local and a bi-local term. The former can be written as in terms of the stress tensor On the other hand, the bi-local term couples the field between different intervals.
Let us now proceed to the case of a chiral fermion on the torus. As is customary, we take the periods to be 1, τ with (τ ) > 0, such that the nome q := e iπτ is inside the unit disk. We restrict to purely imaginary modulus τ = iβ, where β is the inverse temperature -the general case can be recovered by analytic continuation. For simplicity, we move to radial coordinates w = e iπz .
Since we are dealing with fermions, the correlator G(u, v) with u = e iπx and v = e iπy is either periodic (Ramond; R) or anti-periodic (Neveu-Schwarz; NS) with respect to either of the two periods of the torus. We shall restrict to the "thermal" case, with NS periodicity with respect to τ . Combining this with the requirement to reproduce the UV correlator G UV (x, y) = [2πi(x − y)] −1 on small scales, this fully determines the standard Green's functions [41] where η(q) and ϑ ν (z|q) are the Dedekind eta and Jacobi theta functions (see supplemental material).
Here, the superscript ν = 2, 3 = (R, NS), (NS, NS) labels the different spin-structures, and we introduced a regulator in order to treat the distribution G ν as a function. The sign of depends on the chirality-without loss of generality, we choose > 0.
With the notation settled, we now go back to the integral equation (6). In radial coordinates, it reads with A := e iπV being the entangling region. The key observation of this paper is that (9) resembles the result of a contour integral, involving simple poles and branch cuts. Thus the strategy to solving (9) is to recast it as a contour integral.
To this end, we start by listing a set of sufficient properties that F ν V must possess in order to solve this equation: (A) Periodicities. First, it must have the same periodicities in the w argument as G ν , such that G ν F ν V is well defined on the torus. The reason is that doubly periodic functions have vanishing residue along the boundary γ of any fundamental region (see Fig. 1): Our aim is now to rewrite this in the form of (9).
(B) Location of poles and branch cuts. The next property we demand is that together with a branch cut along the entangling region A, which we specify below. Everywhere else it must be analytic. Note that, similarly to G ν , we need to introduce a regulator > 0 for the pole of F ν V .
FIG. 1: The complex plane analysis in the argument. The black solid line is the entangling region-here for simplicity two intervals. The blue line represents the contour of integration γ in (10), which leads to the residues evaluated along the green dotted curves.
If these conditions are met, a simple residue analysis shows that (10) reduces to where we made the regulators explicit and A denotes a snug path around the cut on A as depicted in Fig. 1.
(C) Residues. This last integral decomposes into three contributions: one along A just inside the unit circle, one along A just outside the unit circle, and contributions from the boundary points α n = e iπan , β n = e iπbn of A as can be seen from Fig. 1. Our final requirements on F ν V are that the residues at ∂A vanish, while F ν V has to have a multiplicative branch cut along A: at every point, the ratio of the function just above and below the cut is a fixed number The solution to (12) in the plane is familiar: F (z) = z m with m ∈ Z possesses such a cut. Below we find the analogue of this on the torus.
If properties (A),(B),(C) are satisfied, it is easy to show that such an F V indeed solves the problem: our contour equation (11) becomes exactly the original singular integral equation (9). The requirement that the residues on ∂A vanish is equivalent to demanding that the modular flow behaves like Rindler space in the vicinity of ∂A. This is analogous to the derivation of the black hole temperature by the smoothness condition at the horizon.
In the supplemental material, we explicitly derive F ν V satisfying all of the above assumptions. The general procedure is as follows: 1. Start with the standard solution for the requirement of a multiplicative branch cut (12) on the cylinder [36].
2. Average over all fundamental domains in the direction of τ . This yields a quasiperiodic function.
3. Multiply with a slightly modified form of the Green's function (8) to turn the quasiperiodicity into a periodicity and introduce the correct pole.
We are now in position to state one of the main results of this paper: the resolvent for a finite union of disjoint intervals on the torus, V = ∪ N n=1 (a n , b n ). The exact expression lives in the complex plane, but is vastly simplified along A. Introducing the shorthand notation where L is the total length of V , our result is with h defined in (12), and Some comments are in order. The term in the second line of (14) is the complex power of a quotient, which introduces the required branch cut along A. This function is quasi-periodic, acquiring a factor of λ 2 when translated into the next fundamental domain. The first factor resembles the propagator (8) and introduces the desired pole, as described above. Additionally, the extra factor of λ in the argument of ϑ ν is there to precisely cancel the quasi-periodicity of the second term. This allows the product G ν F ν V to be exactly doubly periodic, as required.

MODULAR HAMILTONIAN
Finally, now that we have found the resolvent R ν V , we can go back to (3) to obtain the modular Hamiltonian For the sake of readability, we shall keep P implicit for the rest of this paper. Equation (16) implies that the δ-terms from (5) drop out in (3), yielding The main characteristic of (17) is that the integrand is highly oscillatory and divergent around ξ = 1/2. Indeed, notice that when ξ → 1/2 the prefactor in (17) diverges quadratically while F (ξ) vanishes linearly but oscillates wildly due to the last factor in (14). However, this behaviour is well understood in the theory of distributions, and in this sense the expression (17) is well defined and closely related to the Dirac delta.
In the supplemental material, we evaluate (17) analytically. Here we will simply quote the result, but the main steps in the derivation are the following: 1. Change variables to isolates all the infinite poles along the negative axis, which then lie in successive fundamental domains. (17) by placing a contour that includes increasingly many poles, and express it by residues.

Regularize
3. Use the quasiperiodicites of ϑ ν to bring every pole to the fundamental region, expressing (17) as a highly oscillatory function with a divergent prefactor.
4. Remove the regulator, leading to standard Dirichlet kernel representations of the periodic/antiperiodic Dirac delta.
The final expression for the modular Hamiltonian depends on the spin sector. Let us focus on the results for a single interval. Both sectors ν = 2, 3 have a local and a bi-local term. The local term is identical in both cases and takes the form with the entanglement temperature where Ω V is as defined in (15) and the function f (x) is fixed by requiring that K loc is hermitian. Note that the expression (18) is equivalent to the more familiar Rindlerlike representation (7). The bi-local term represents the central result of this paper and shows a novel feature: In both sectors, it involves a coupling between an infinite but discrete set of points, and is given by (20) where the sign ± corresponds to ν = 2 3 . Here, we used the function µ(x, y) = 1 2πL log which will play an important role in the analysis below.
Note that K ± bi-loc couples pairs (x, y) which are solutions of Since µ(x, y) is monotonic in y and diverges at the endpoints, eq. (22) possesses a unique solution for every k, as shown in Fig. 2. Solutions accumulate near the endpoints. In the next section, we analyse the above expressions and discuss their physical meaning. A summary of the results is presented in table I.

DISCUSSION
In this paper we computed the modular Hamiltonian of chiral fermions in a thermal state on the circle, reduced to an arbitrary set of disjoint intervals.
Our main result is that for arbitrary temperature, the modular Hamiltonian contains a local term, as well as an infinite number of bi-local contributions, even for a single interval. Let us now analyse the bi-local terms in more detail. Inserting the kernel (20) back into (2), the bi-local modular Hamiltonian reads As depicted in Fig. 2, the x k (x) are an infinite set of points within the interval, solutions to equation (22). The bi-local coupling α(x, x k ) has dimensions of energy and is given by Although determining the exact location of the x k is difficult, two properties are simple to extract: • The infinite set of x k accumulate near the endpoints of the interval. Indeed, since µ diverges there, there is an infinite number of solutions near the boundaries, located at and similarly near b.
• Their contributions vanish as they approach the endpoints. Using (20), the coupling in (23) goes as The energy scale of α(x, x k ) is set by the temperature β −1 , whereas the fall-off is determined by the length of the interval L. Interestingly, the strength of the nonlocal couplings appears to be "redshifted" due to their proximity to the local Rindler horizons located at the endpoints.
As a next step, let us see how to recover the known results at very high [34] and low [36] temperatures. We start with the high temperature limit β → 0. One easily sees from (19) that the local term goes as the inverse temperature, β(x) ∼ β, as expected. On the other hand, as depicted in Fig. 2, the bi-local contributions (20) all approach the endpoints, where they vanish exponentially.
Moving now to the low temperature limit β → ∞, the entanglement temperature (19) approaches the well known result for the cylinder [36] lim β→∞ β(x) = 2π The bi-local contributions however behave remarkably. As can be understood from Fig. 2, as we lower the temperature, the curve gets increasingly steep. Thus, the solutions to (22) form a partition of the interval which becomes denser and denser in the limit β → ∞. Now, recall that the modular Hamiltonian must always be thought of as a distribution, i.e. as integrated against regular test functions. In this limiting procedure, the solutions to (22) "condense" in the interval, and it can be shown that the sequence of Dirac deltas in (20) reproduce precisely the definition of a Riemann integral. Indeed, one can show that in this sense (20) becomes completely nonlocal in agreement with [36], whereas lim β→∞ K − bi-loc = 0 due to the oscillating (−1) k .
The previous analysis provides a new insight into the structure of fermionic entanglement: At any finite temperature, non-locality couples a given point only to an infinite but discrete set of other points. The characteristic scale needed to resolve this discreteness goes as 1/β. Hence, continuous non-locality emerges strictly in the limit of zero temperature. We summarize the structure of the modular Hamiltonian in table I.  (18)- (20). The local and non-local terms at low temperature (β → ∞) are given in (26) and (27).
For multiple intervals, the only difference is that (22) now possesses one solution per interval for a given k, including the non-trivial (x = y) solutions for k = 0. In the low temperature limit, these extra terms yield precisely the bi-local terms of [35,36].
During the final stage of this project, related results were independently reported in [42] and [43]. Eqs. (145) and (146) of [42] give the modular flow of the correlator. The generator of this flow corresponds to the expectation value of our result for the modular Hamiltonian. Finally, the versatility of the resolvent method has allowed to compute the associated entanglement entropy [44], and can also be used to study other quantities related to the entanglement spectrum.

ACKNOWLEDGMENTS
We are very grateful to D. Blanco and G. Pérez-Nadal for collaboration in the initial stages of this project.
We thank J. Camps, B. Czech, M. Heller, H. Hinrichsen, C. Northe and G. Wong for helpful discussions on the subject. PF is financially supported by the DFG project DFG HI 744/9-1. The Gravity, Quantum Fields and Information group at AEI is generously supported by the Alexander von Humboldt Foundation and the Federal Ministry for Education and Research through the Sofja Kovalevskaja Award. IR also acknowledges the hospitality of Perimeter Institute, where part of this work was done.

Conventions
We work with the Dedekind eta and Jacobi theta functions defined as: The resolvent from the propagator Let us give a quick derivation of (3). The basic idea is that given a holomorphic function f (z) and an operator G, we can determine f (G) by using Cauchy's integral formula for each of the eigenvalues of G.
In our case, we wish to find K = − log[G −1 − 1]. Since the spectrum of K is real, this equation tells us that the spectrum of the propagator G is a subset of the open interval (0, 1). Consider now a specific eigenvalue g of G.
We can use Cauchy's integral formula to find where γ g is a contour that encircles both g and 1−g. Notice that log z possesses a branch cut, which we choose to place along the negative real axis -the contour γ g must not intersect this cut. Then, since the function is holomorphic everywhere else in the plane, we can freely deform γ g such that integration only has to be done along the branch cut and a circle at infinity. The latter contribution however vanishes since the integrand is bounded by z −2 log z for large z, and we obtain with log ± referring to the values just above and below the branch cut. After evaluating the last bracket to 2πi, we can change variables to ξ = 1/2 − z to find As this formula is valid for every eigenvalue g of G, it also holds as an operator statement, hence we arrive at (3).

Deriving the resolvent
In this section, we derive the solution F ν V (given in (14)) to the singular integral equation (9). Let us start with the functions [36] ω n (w) := sin(π(a n − z)) which provide the correct branch-cut on the cylinder. Choosing the branch cut of the logarithm along the negative real line, we see that for w ∈ A. Note that (28) is not well defined on the torus since it transforms non-trivially under w → qw. We shall remedy this by defining log Ω n (w) := k∈Z log(ω n (q k w)) − log(ω n (q k )) , (30) where the second term in the brackets is to ensure absolute convergence. We made the logarithms explicit in order not to break the behaviour (29) at the branch cut. At first sight, Ω(w) seems doubly-periodic: by construction, ω is periodic with respect to the spatial circle, and now we sum over all translations along imaginary time. However, Ω(w) has a non-vanishing residue within each fundamental region due to the branch-cut, and thus cannot be elliptic. Instead, it turns out to be quasiperiodic, as is seen by putting a cutoff in the sum (30), and then computing Ω(qw). Then, the series acquires a prefactor originating in lim k→±∞ ω n (q k w) = β n α n
This yields the quasi-periodicity where λ is defined in (13). To cancel off the quasiperiodicity and to introduce the desired pole, we multiply with a combination of theta functions. We find where we made our choice of branches in Ω ih n explicit (our choice is such that the residue evaluation in the main body of the paper does not cross the branch cut). Eq. (31) now solves (9).
Finally, let us rewrite this in terms of more familiar elliptic functions. Note that (29) and, now that the numerator and denominator are on the same side of the branch cut, we can move the product into the complex power to find After some algebra and an application of the Jacobi triple product [45], this simplifies the solution (31) to (14) with Ω V from (15).

Deriving the modular Hamiltonian
In this section, we provide the mains steps to evaluate the integral expression (17) for the modular Hamiltonians. We restrict to purely imaginary τ = iβ-the general case can be restored by analytic continuation. Let us first change the variable of integration from ξ to Λ := λ 2 = e 2πLh = ξ + 1/2 ξ − 1/2 L , such that (17) turns into where we use the shorthand notation µ := 1 2πL log To evaluate the above integral, note the following two facts: • Since we merged the two occurences of F ν V (±ξ) in (17) into a single integral, integration has to be done symmetrically with respect to ξ ↔ −ξ, i.e., Λ ↔ Λ −1 .
This requires that we introduce a symmetric regulator r (Λ) = r (Λ −1 ) to tame the integral, allowing us to evaluate it via standard complex analysis methods. We choose r (Λ) := (1 + ) 2 (Λ + )(Λ −1 + ) (32) to otain The integral (33) can now be evaluated using contour integration. To this end, consider the integral where the contour γ is as depicted in Fig. 3. The circular contributions vanish due to the falloff of the regulator. Choosing the branch cut of Λ iµ along the positive real axis, we see that two remaining horizontal contributions yields two almost identical terms, differing differ only by a global prefactor of −e −2πµ . We thus find By Cauchy's theorem, I ν can also be expressed as a sum over residues, yielding a series expression for K ν V . We shall do this explicitly for ν = 3 and briefly mention the differences for ν = 2, 4 at the end.
The poles of the integrand are of two types (see Fig. 3): Two come from the regulator (32), located at Λ → − and Λ → − −1 . The other (infinitely many) poles come from the poles of the 'propagator-like' term. As can be seen from either the Laurent expansion of this term (see section below) or directly from the Jacobi triple product, we have the leading divergences at Λ → −q 2k+1 for k ∈ Z. Keeping in mind that the negative sign of poles always has to be written as e +iπ due to our choice of branch cut, this yields Let us have a look at the series in the last line: Using the Laurent expansions below, this can be rewritten as We choose the cutoff to be = q 2m with very large m ∈ Z to avoid the poles at q 2k+1 , so that we only deal with simple poles. Then, putting everything together into (36) and using the quasiperiodicities of ϑ 4 , one finds As already stated above, this limit must be understood in the sense of distributions. We see that K 3 V contains essentially two factors: the term involving the sine function is highly oscillatory for m → ∞, except at solutions of As a distribution, it vanishes when integrated against any regular test function. However, the remaining factor P (x, y) is not regular since it has poles, and thus we must examine its behaviour in their vicinity, which will lead to finite contributions. These poles coincide precisely with the solutions to (37), which are of two kinds: the trivial solution x = y will lead to a local term, while the other solutions x = y will give bi-local contributions. Let us start with the latter. Close to these solutions, a straightforward calculation shows that P (x, y) ∼ iπ L sinh πµ(x, y) 1 sin(π[x − y + βµ(x, y)]) .
Now we turn to the solution x = y, which is special as it leads to a second order pole in P . Similarly to before, in the vicinity of that solution, P (x, y) takes the form which together with the oscillatory term leads to − iβ L δ(x − y + βµ(x, y)) x − y .
Note that we did not need to consider the terms with k = 0 as in (38) since we only deal with the solution x = y. As a last step, we use the methods from [35] to rewrite the singular fraction (40) as β L [i∂ x + f (x)]δ(x − y) 1 + β(∂ x µ)(y, y) , where f (x) is fixed by hermiticity. Now we focus on the case of a single interval. Again we begin by considering on the bi-local terms. Since µ(x, y) is monotonically increasing with respect to x in the interval, eq. (37) has a unique solution for each k ∈ Z. In particular, note that for k = 0, the solution is x = y. Since we already consider this contribution separately in (41), we can explicitly exclude it from the series (39). The final expression for the modular Hamiltonian for a single interval is then given by the sum of (18) and (20). Finally, replacing (20) in (2), the bi-local modular Hamiltonian takes the form (23).
An analogous calculation holds for ν = 2, with one small adjustment: Since the poles of the Laurent expansion are instead located at −q 2k , we obtain the periodic version of the Dirichlet kernel in (38). The rest of the calculation is identical.

LAURENT EXPANSION
To better understand the location and behaviour of the poles of the "propagator-like" terms in (14), we derived their Laurent expansions. The coefficients may be computed as contour integrals which vastly simplify due to the quasi-periodicities of the theta functions. In the fundamental domain |q| 1/2 < |w| < |q| −1/2 , the result then takes the form of Lambert series