Entanglement, Replicas, and Thetas

We compute the single-interval Renyi entropy (replica partition function) for free fermions in 1+1d at finite temperature and finite spatial size by two methods: (i) using the higher-genus partition function on the replica Riemann surface, and (ii) using twist operators on the torus. We compare the two answers for a restricted set of spin structures, leading to a non-trivial proposed equivalence between higher-genus Siegel $\Theta$-functions and Jacobi $\theta$-functions. We exhibit this proposal and provide substantial evidence for it. The resulting expressions can be elegantly written in terms of Jacobi forms. Thereafter we argue that the correct Renyi entropy for modular-invariant free-fermion theories, such as the Ising model and the Dirac CFT, is given by the higher-genus computation summed over all spin structures. The result satisfies the physical checks of modular covariance, the thermal entropy relation, and Bose-Fermi equivalence.


Introduction
The study of entanglement measures in quantum field theories has provided a rich set of results that illuminate these theories as well as their holographic duals (for a comprehensive review, see [1]). In this work we focus on some of the simplest known quantum field theories, namely free conformal field theories in 1 + 1 dimensions, and address some puzzles that arises when one computes Rényi and entanglement entropies for systems of finite size and at finite temperature.
The first such computation was performed in [2] in the case of (1 + 1)-dimensional fermions with a single entangling interval. This calculation was performed using the replica trick [3][4][5] in a standard way and computing two-point functions of twist operators on the torus. For fixed fermion boundary conditions (spin structure) on the torus, one finds a result in terms of product of Jacobi theta functions. Let us take the size of the spatial circle to be L and the entangling interval to be , and define a scaled interval z 12 = L that lies between 0 and 1. Then the entaglement entropy S EE (z 12 ) satisfies the relation [2]: (1.1) This "thermal entropy relation" was originally proposed in [2] from a holographic point of view and later derived directly in CFT in [6], see also [7][8][9]. The free-fermion calculation was subsequently streamlined and generalised [10] to the massive case, still for a fixed spin structure.
Path-integral computations of the Rényi entropy at finite size and temperature give an answer that is (at least locally) analytic in the modular parameter τ of the torus and the interval size z 12 . Thus it is natural to ask if the result is modular covariant 1 . In this context the following puzzle was raised in [12]. Firstly it was observed that, although the replica partition function for any fixed spin structure is not modular covariant with respect to the full modular group of the torus, one can obtain a modular covariant answer by summing over all four torus spin-structures for the fermions. While this sum over torus spin-structures may seem to be a natural observable from the torus point of view, it does not satisfy Bose-Fermi equivalence, which equates the free fermion theory (after summing over spin structures) to a suitable free boson theory. Furthermore, this summed result does not obey the thermal entropy relation above.
In this paper we resolve these puzzles. Two concrete examples, on which we will focus for the most part in this paper, are the Ising model i.e. one Majorana fermion summed over spin structures, and the Dirac CFT i.e. one Dirac fermion summed over spin structures.
Each of these is a well-defined 2d CFT with a known partition function, and more general examples can be found in [13]. The Rényi/entanglement entropies of these theories exhibit 1 Physically the length of the entangling interval is real, and the modular transformation S transforms z12 from real to imaginary values. Thus the modular transformation of entanglement entropy is a "temporal" version of it. However, the path integral calculation is naturally performed for a generic complex interval on the 2d torus spanned by ordinary space and Euclidean time. The result is holomorphically factorised in the interval size z12(z12), and it is meaningful to ask in this analytically continued setting if it is modularcovariant. Note that in a different context, that of multi-interval entanglement at zero temperature, it has also been emphasised in [11] that entanglement entropy should satisfy both modular invariance and Bose-Fermi duality.
the above puzzles very precisely. For example the Dirac CFT is dual to a free boson theory with R = 1. However, the free Dirac fermion with any fixed spin structure is not equivalent to a free boson theory. The thermal entropy puzzle is also illustrated clearly in this system. The twist-operator computation exhibits an ambiguity in whether the spin structure should be summed over before or after taking the product of θ-functions across replicas. As shown in [12], neither way of performing the sum satisfies the thermal entropy relation. A strong form of the thermal entropy relation, as explained in [6,8], says that the n th replica partition function Z n (z 12 , τ ) should reduce in the limit z 12 ∼ 0 to Z 1 (τ ) n (up to the prefactor that encodes the singularity induced by the collision of the two twist operators), while in the limit z 12 ∼ 1 it should go over to Z 1 (nτ ). What was shown in [12] is that the former prescription only agrees with this prediction at small intervals z 12 ∼ 0, while the latter prescription only agrees with it at large intervals z 12 ∼ 1. This means that the twist-operator computation, at least in the form currently known, is inadequate to compute the (finite size, finite temperature) Rényi entropy of the free Majorana (Ising)/free Dirac CFT. Subsequent to [12], issues regarding summing over fermion spin structures when computing Rényi entropies have been discussed in different contexts in [14][15][16][17][18].
The replica partition function has also been computed at finite size and temperature for free bosons compactified at an arbitrary radius R in the target space. The correct answer appears in [19] (building on previous work in [20], in turn based on the orbifold twist-operator computations of [21]). In [19] the answer was shown to satisfy the thermal entropy relation, while in [12] it was shown to also be modular covariant. This result is essentially the higher-genus partition function of the free boson theory on a Riemann surface of genus n, where n is the number of replicas. This Riemann surface has a fixed and rather special period matrix Ω determined by two complex variables that are the physical parameters of the original problem: the modular parameter τ of the original unreplicated torus and the length z 12 of the sub-system.
We propose here that the only correct way to compute the Rényi entropy of modularinvariant free fermionic theories at finite temperature and size is through the higher-genus relica partition function approach 2 . The full partition function is computed on the genus-n replica surface as a sum over all the 2 2n spin structures 3 of this surface, α = (α 1 , α 2 , · · · , α n ) and β = (β 1 , β 2 , · · · , β n ) whose entries are independently chosen to be any integer or halfinteger. We show that the final answer satisfies (i) modular covariance on the original torus and (ii) the thermal entropy relation. Bose-Fermi equivalence is relatively trivial because 2 Multi-interval Rényi entropies for free fermions at zero temperature have been studied in this way in [22,23] and WZW models have been studied in [24,25]. 3 The partition function vanishes for odd spin structures, so only the 2 n−1 (2 n + 1) even spin structures contribute to this sum.
the higher-genus answer is already known to satisfy it for any period matrix Ω, and here we have simply specialised Ω to the replica surface. Furthermore, if we focus on special spin structures of the higher genus partition function, the resulting expressions turn out to have interesting symmetry properties, in that they transform as Jacobi forms in the variables (τ, z 12 ).
Our proposal leads to some non-trivial checks. As we have indicated above, and will elaborate in what follows, for certain fixed spin structures the n-th Rényi entropy can be calculated in two distinct ways-as a partition function on the genus-n replica surface, and as a two-point function of the n th twist operator on the original torus. These spin structures are the special "diagonal" higher-genus spin structures that obey the replica symmetry, namely those of the type α diag = (α, α, · · · , α), and β diag = (β, β, · · · , β), where (α, β) is the spin structure on the original torus 4 . The higher-genus answer is expressed in terms of Siegel Θ-functions with the special characteristics α diag β diag . The twist-operator approach involves computing two-point functions on the original torus, and the answer is expressed in terms of Jacobi θ-functions with spin structure (α, β). We find that these two calculations are equivalent only if a very precise identity is satisfied, which relates Siegel Θ-constants evaluated on the special replica period matrix to a product of Jacobi θ-functions of the variables (τ, z 12 ).
We study this proposed identity in detail, showing independently that each side has the same periodicity in the entangling interval z 12 as well as the same transformations under modular transformations of the original torus. As a result both sides transform as Jacobi forms of the same weight and index. The identity itself turns out to be quite non-trivial, and we check it by expanding each side in powers of z 12 . At low orders we demonstrate the identity explicitly for all n. For higher orders, we consider the special case of two replicas n = 2, and evaluate the difference of the two sides up to 40th order in z 12 to find an equality result in each order. This constitutes strong evidence that the higher-genus and twist-field computations agree.
Next we elaborate on a certain interesting aspect of the results. As defined, the scaled interval length z 12 takes values between 0 and 1. Thus, in principle the cut should lie within a fundamental region of the original torus. Once we analytically continue in the cut-length z 12 , however, there are other paths connecting the same two endpoints that wind around the two cycles of the torus. The inequivalent paths turn out to be those that wind 0, 1, · · · , n − 1 times around either cycle. The fundamental region of our problem is therefore a torus of sides (n, nτ ). We find, in accordance with this intuition, that the partition function has a periodicity 5 under shifts of z 12 by the lattice C/(nτ Z + nZ).
The plan of the paper is as follows. In Section 2 we discuss the special higher-genus surfaces arising from the replica trick and study their partition functions. We then present our proposal of equality of the two ways of computing these partition functions-the direct higher-genus method and the twist-operator method. In Section 3 we discuss various checks of our proposal for arbitrary n, including the symmetry properties and the small-interval expansion. In Section 4 we focus on the n = 2 case and verify our proposed relation to high order in an expansion in z 12 . In Section 5 we sum over spin structures and show that the answer obeys the various physical properties that we expect. In Section 6 we summarize our results and conclude. In the two appendices we discuss, respectively, the periodicity and modular properties of the higher-genus and twist-operator expressions.
2 Higher-genus surfaces and replica partition functions

Period matrix of the replica higher-genus surface
Let us start by summarising the methodology of computation of the replica partition function for a CFT at finite size and temperature, i.e. on the torus. As is well-known, in such a case the replica trick gives rise to a Riemann surface made up by joining n copies of the torus (where n is the number of replicas) sequentially along the entangling interval. The result is a genus-n surface, though a very special one. By definition, the replica partition function should be simply the partition function of the given CFT on this genus-n surface. On the other hand, one can introduce "twist operators" whose non-local OPE's with the free fermions have the effect of "transporting" them from one replica to the next. The replica partition function is then identified with the correlation function of these twist operators. Next, a diagonalisation procedure, described for example in [5], reduces the problem to the correlation function of twist fields on a single replica surface. In this approach the computation is performed on a torus rather than on the replica surface.
The question we wish to first address is whether the two computations, which we call the "higher-genus calculation" and the "twist-field calculation", are equivalent at least for a fixed free-fermion spin structure. Since there are more spin structures in genus-n than on the torus, this question is not well-defined. Nevertheless, one can ask if the highergenus calculation for the special class of spin structures of the form α diag = (α, α, · · · , α) (and similarly for β) gives the same answer as the twist-operator calculation for a given 5 There is a small subtlety that for even n each fixed-spin-structure partition function transforms into itself under shifts of 2n and 2nτ , but in fact it transforms into a different spin-structure under shifts of n and nτ so that the full summed partition function transforms correctly under shifts of the lattice C/(nτ Z + nZ).
spin structure (α, β) on the torus. Note that the two calculations are quite distinct. The higher-genus one uses general results about free-fermion partition functions on higher-genus surfaces [26,27]. One takes the appropriate result and inserts the period matrix of the replica Riemann surface, which as we will soon see is quite special. On the other hand the twist-operator computation uses the torus correlation function of twist operators. With the higher-genus spin structures restricted as above, we will find a precise equivalence between them. This result follows from a nontrivial identity between a genus-n Siegel Θ-constant evaluated for a specific subclass of period matrices, and a product of genus-1 Jacobi θfunctions. We will show that both sides of the identity are in fact Jacobi forms [28] on an n-fold/2n-fold cover of the original torus (depending on whether n is odd or even).
Riemann surfaces of genus n are parametrised by their period matrix Ω which is in general a complex symmetric n×n matrix. As is well-known, this has more parameters than necessary to describe a general Riemann surface of this genus. The problem of determining which subset of Ω's correspond to a Riemann surface is called the Schottky problem.
However the situation of interest to us is much simpler, as we are only interested in the special Riemann surfaces that arise via the replica trick. The surface defined by taking n replicas of a torus with a single interval has only two independent moduli, z and τ , where τ = i β L is the modular parameter of the original torus and, as already indicated, z 12 = L is the relative length of the entangling interval. Thus, we would like to express Ω ij as a function of z 12 and τ .
The desired answer follows from the constructions of cut differentials for the Riemann surface of interest [19][20][21]. These are given, for example, in the Appendix of [19] 6 : where n is the number of replicas which equals to the genus of Riemann surface and k = 0, 1, 2, · · · , n − 1 labels the linear independent differential cuts. Selecting a basis of cycles A a , B a , we have: The period matrix of the Riemann surface is then Ω = B · A −1 .
So far we have not specified a basis of cycles. However there is something special about the genus-n Riemann surface obtained by gluing n tori sequentially along a cut, namely the fact that the gluing procedure does not introduce any new handles to the surface. It simply 6 The arguments of our cut differentials are shifted with respect to those of [19] so that ω here (z) = ω there z + (1 − k n )z1 + k n z2 . This does not, of course, affect the periods.
connects n genus-1 surfaces into a single genus-n surface. This is in contrast to the surfaces obtained by gluing several complex planes along a pair of cuts [29], relevant to the case where the entangling region is made of two disjoint components (at zero temperature and infinite spatial size). In the latter case it is in fact the gluing procedure that introduces the handles to the resulting surface. Therefore in that case, any basis of cycles on the replica Riemann surface must necessarily involve the cuts. For our case, each component torus of the replica surface already has a pair of canonical cycles, the usual (A, B) pair for a torus.
These continue to be valid cycles of the glued surface of genus n, and it turns out very convenient to choose them as a basis for the latter. Accordingly, from now on A a , B a will be taken to be the cycles of the a-th torus component of the glued replica Riemann surface.
This choice will considerably simplify the analysis of the problem.
By setting k = 0, we see that A 00 = 1, B 00 = τ . We can now relate all the entries A ak to A 0k , and likewise B ak to B 0k . For this, notice that ω k picks up a phase α k where α = e 2πi n , when we go from one replica to the next one above it. It follows that: where: We can think of the above as products of the matrix M with entries M ak = α ak with Then we have: where we have defined the diagonal matrix C = diag(C k ) with components 7 : From the discussion above, we know that C 0 = τ . The inverse of M (α) is 1 n M (α −1 ) and one has: (2.8) 7 In the notation of [19], this corresponds to C k = Using the property that C k = C n−k , which is easily verified, we can check that the above matrix is symmetric as it should be. Hence it can equivalently be written as: As examples, for n = 2, 3 we have: (2.10) For every n the period matrix satisfies the identity that the sum of all elements in a row For future use, we note here some additional identities that it satisfies: (2.12)

Relationship between higher-genus partition function and twist-operator computation
We would now like to describe the relationship between the two ways of computing the replica partition function: one as a higher-genus partition function evaluated on the restricted family of period matrices defined above, and the other as a correlation function of twist operators on the torus.
For the former, we first write the genus-n Siegel Θ-function: where the characteristics α = (α 1 , α 2 , · · · , α n ) and β = (β 1 , β 2 , · · · , β n ) are two n-component vectors whose entries are independently chosen to be any integer or half-integer. The independent choices are 0, 1 2 , and all other choices can be related to these. This theta-function is one of the factors in the free-fermion partition function in higher-genus. In the context of free fermions, the characteristics arise as boundary conditions along different cycles of the Riemann surface, namely spin structures.
The genus-n fermion partition function for a single Dirac fermion with arbitrary highergenus spin structure α β is [26,27]: while for a single Majorana fermion, it is: where C is a spin-structure-independent factor related to the determinant of an (anti-) holomorphic differential operator on the surface (see Eq.(5.13) of [26]) 8

. The full Dirac or
Majorana partition function is a sum over all 2 2n spin structures arising on the higher-genus Riemann surface: We see that besides the Θ-function, the replica partition function has an additional factor of a power of |C|. To compute it for our case, we use the fact that this determinant is independent of spin structures, and is the same one that appears in the bosonic partition function. Indeed, it was proved long ago [26,27] that Bose-Fermi equivalence holds on arbitrary Riemann surfaces. Now the free boson replica partition function [19,20] is expressed as a higher-genus boson partition function evaluated on the period matrix Ω of the previous section. It depends on the radius R at which the boson is compactified. For R = 1 one can perform standard manipulations to reduce the free boson partition function to a sum over Θ-functions with characteristics α, β times some other terms. The result is: (2.17) This bosonic partition function can be understood as arising from a classical part and a quantum part. The classical part is simply the higher-genus theta function, so the quantum 8 As pointed out in [26] this quantity depends not only on the moduli but on the metric of the surface because of conformal and diffeomorphism anomalies. For us the precise metric is determined by the replica construction to be flat everywhere except at the end points of the cuts.
part can be identified with the factor 1 2 g |C| 2 . We can thus read off this factor to be: For a Majorana fermion one simply takes the square root of the above expression.
The question is now whether the higher-genus approach leading to Equation (2.14), (2.15) and the torus twist-operator approach leading to Equation (2.19) are equivalent. As it stands, however, this is not a well-posed question because the former depends on 2 2n spin structures while the latter has just 4. Thus we restrict the higher-genus expression to a fixed spin structure α diag = (α, α, · · · , α) and β diag = (β, β, · · · , β). These are very special spin structures that are taken to be the same over each "handle" of the genus-g surface.
They are certainly not the only ones that contribute to the full modular-invariant higher- Both the objects in question are squares of locally analytic expressions. Thus we take the holomorphic square root χ, so that Z = |χ| 2 , and examine whether the two expressions with Ω defined in Equation (2.9), and are equal to each other.
Cancelling out common factors, we can restate the above considerations in terms of the higher-genus expression 3 Some checks of the higher-genus partition function/twist-operator equivalence In what follows, we provide evidence for Equation (2.24) which is a nontrivial equality between two well-defined mathematical functions. To start with, we plot them numerically for n = 2. We fix the modular parameter to be τ = i and the spin structure to be (α, β) = (0, 0), and let z take real values from 0 to 1 2 . The plots are shown in Figure 1. The proposed equality is exact, as our subsequent mathematical analysis will show, hence the slight deviations visible near z = 1 are attributable to an inaccuracy of the numerical plot rather than of the equality itself. This in turn arises from the fact that the cut is 0 < 2πz < π (i.e. 0 < z < 1 2 ). This is for Dirac fermions (for Majorana fermions one would take the square root of both sides). differential involves a square root and one has to choose the integration contour to avoid the branch cut. The numerical programme becomes inaccurate as the contour approaches one edge of the cut.
Our main focus in the rest of this paper is to provide mathematical evidence for the equality using known properties of the functions on both sides. There are some basic checks that must be satisfied if it is to hold. A zeroth check is that both sides are manifestly equal for n = 1. Next, when there is no cut, the genus n surface factorizes into n identical genus-1 surfaces. In accordance with this, when z 12 = 0 we have A 0k = 1, B 0k = τ and Ω ab = τ δ ab , so that Now note that if n is odd then χ t is zero for the spin structure because this is an odd spin structure on the torus and the k = 0 term in the product therefore vanishes. At the same time, χ g vanishes because for odd genus this corresponds to an odd spin structure.
However for even n there is no k = 0 term in χ g and it does not have to vanish even for the spin structure . And from the higher-genus point of view, this spin structure is now even and therefore χ t also does not have to vanish. Thus some simple checks are satisfied.

(3.3)
We see that the periodicities (3.3) are not those of the torus but of its n-fold cover. As explained in the introduction, this is due to the fact that there are non-trivial paths for the cuts connecting the two end-points z 1 and z 2 which wind around the two cycles of the torus. A path that winds around n times around either cycle of the original torus is equivalent to a closed cycle on the higher-genus surface, and therefore the inequivalent paths are those that wind 0, 1, · · · , n − 1 times around either cycle. In the context of the Ising model these winding paths can be identified with disorder operators [30].

(3.7)
We see that our quantity F (Z, τ ) satisfies these equations, upto a phase and also up to a change in the spin structure (α, β) 9 . From the transformation properties of Equations (3.5) and (3.6) we deduce that χ g and χ t have weight and index: For even n, the above discussion has to be slightly modified. Under shifts of n and nτ , both χ g and χ t change their spin-structure (in the same way). Therefore to find their quasi-periodicity at fixed spin-structure, we must double the shifts. Indeed we show in Appendix A that, for even n, χ g and χ t are both periodic under z 12 → z 12 + 2n and quasi-periodic with the same prefactor under z 12 → z 12 + 2nτ .
We note here that the function χ t is manifestly holomorphic in z 12 , which allows us to write a double Fourier expansion: From the properties of Jacobi theta functions we also see that n only takes positive values in some one-dimensional lattice depending on n and the spin structure (i.e. n, r need not be whole integers). Combined this observation with the transformation properties above, we see that χ t is really a weak Jacobi form in the sense of [28].

Small-interval expansion
A stronger check is to expand both sides in powers of the interval size z 12 . We have already seen that χ g and χ t agree at z 12 = 0. As we see below, the coefficient of z 2 12 vanishes on both sides, and the first non-zero coefficient is at O(z 4 12 ). In this subsection we consider both n (the number of replicas) and the spin structure to be arbitrary (and correspondingly we suppress the spin structure label). We compare the expressions χ g and χ t to order O(z 4 12 ) in the interval size, and find non-trivial agreement for all n.
On expanding χ t from Equation (2.23), we get: where denotes the derivative with respect to the first argument. Here and in the following, for simplicity of notation, whenever the first argument of the θ function is suppressed it is understood to be 0. Since each Jacobi θ-function is either even or odd as a function of its first argument, alternate terms in the above expansion vanish-though we will carry along all terms in the interest of a uniform notation. For the odd spin structure there are also vanishing θ(0|τ ) factors in denominators and in front of the full expression, of course the two cancel each other out. In this way we can use the same formulae for all spin structures.
Let us now switch to the variable Z = z 12 /n defined earlier. Working to quadratic order in Z, the last factor above is: (3.11) Note that the term of order Z vanishes.
Thus we have: Now we would like to compare χ t with χ g defined in Equation (2.22) for general n and general spin structures, up to O(Z 2 ). For this, let us recall the cut differential and express it in terms of Z: where the τ -dependence has been suppressed to simplify the notation. Notice that this is invariant under the simultaneous transformation k → n − k and Z → −Z. We expand this to second order in Z and find: (3.14) Again, the term of order Z vanishes. The next step is to compute the integrals: Here we used the identities: (3.16) It follows that: Next we compute the matrix Ω from this using Equation (2.8): (3.18) Defining: we may write: Now we are in a position to evaluate the function χ g in Equation (2.22) to second order in z 12 . To this order, we have seen above that A 0k = 0 for all k. Thus to this order, we have: where we recall that α diag := (α, α, · · · , α) and similarly for β.
Expanding the Θ function we get: The correction term can be written: from which one finds:

(3.24)
It is easily shown that: where the second equation follows from a f (a − b) = 0 combined with the first equation.
Inserting these into Equation (3.24), we find the correction factor to be: 26) in perfect agreement with Equation (3.11).

Verifying the equivalence for n = 2
In this section we focus on the case n = 2 and aim to establish an equality χ t = χ g in a power-series expansion in Z = z 12 /2. We do this by expanding both quantities as a power series in Z, and expressing each coefficient as a function of τ only. The coefficients turn out to be functions of the Jacobi theta constants θ α β (0|τ ), their derivatives, and the Eisenstein series: (4.1) The twist field expression χ t can be written for even spin structures as the product of two theta functions: The odd spin structure We now write this as a power series (suppressing α, β on both sides, and denoting by f (n) (z) the n th derivative with respect to z): Now we us consider the higher-genus expression χ g for n = 2. In this case there are two cut differentials ω 0 = 1 and ω 1 given by Equation (2.1). Since there is only one non-trivial cut differential, we denote it by (with θ 1 the odd Jacobi theta function): The period matrix (2.9) is given by: where C = B/A. Defining B and C via the relations we write the period matrix as: The higher-genus expression is given by: which we want to expand in powers of Z. The expansion only contains even powers of Z as all functions are even functions of Z. First we expand log ω: Expressing the logarithmic derivative of θ 1 in terms of the Weierstrass ℘-function: and using the fact that ω(0, Z) = 1, we obtain the expansion of the cut-differential: Next we define the coefficients of the periods defined in (4.5) and (4.7): In order to compute these coefficients as a power series, we begin by writing the Laurent expansion around 0 of the Weierstrass ℘-function: This expansion implies the following useful equation for every integer j > 0: using which we obtain the periods: Putting this together with the expansion (4.12), we obtain the coefficients A 2n , B 2n .
The first few coefficients in the expansion of A are: and those of B are: (4.18) Finally, expanding C as C(Z) = ∞ 0 C 2n Z 2n , we have the first few coefficients: The Siegel Θ function in Equation (4.9) can be written as: 20) where in the last step we defined u = 1 2 (w 1 + w 2 ), v = 1 2 (w 1 − w 2 ) and used the fact that we are working with even spin structures. To obtain the required power-series, we expand the exponential of the operator above: where we have again dropped the spin structures and arguments of the θ functions to simplify the notation.
The denominator factor in Equation (4.9) can also be expanded: 1 Putting things together, we find: For the O(z 4 ) term to agree with that of Equation (4.3) we need to show that: To do so we start by writing the identity: which follows from the periodicity property of the Weierstrass ℘-function and using its Laurent expansion (4.15) to make a linear combination regular at the origin. (It also follows by differentiating the famous equation (℘ ) 2 = 4℘ 3 − 60G 4 ℘ − 140G 6 .) We remark that similar equations can be found in this manner for higher powers of ℘ and for products of derivatives, for example: (4.26) Now, evaluating Eqns.(4.11) and (4.25) for z successively equal to 0, 1 2 , τ 2 , 1+τ 2 , one gets the following identity valid for all spin structures: Multiplying by 8θ 2 and rearranging the various terms, we obtain Equation (4.24).
tor D z := 1 2i ∂ z , and the Eisenstein series G 2k (τ ), k = 1, 2, · · · , i.e. they are quasi-modular forms on a congruent subgroup of SL 2 (Z) of weight 2k for the expression at O(Z 2k ) (see e.g. [31]). In this paper we do not give a systematic formal proof for the validity of each of these identities, but perform a computational check of these identities. A proof can be constructed by using the fact that the ring of quasi-modular forms is finitely generated.
It is therefore enough to check a finite number of coefficients in the q-expansion in order to prove these identities. The exact number of coefficients depends on the dimension of the space of quasi-modular forms, and a proof can be constructed for each k by using the dimension formula 10 . Using the PARI/GP program [32], we checked that the coefficients of the functions appearing up to O(Z 40 ) each vanish up to O(q 400 ). We consider this convincing evidence that χ t,n=2 and χ g,n=2 are equal at each order in Z.
The physical intuition behind our proposal, as well as our modular forms calculations, suggest that there is a more formal and elegant mathematical proof of these relations. We note that relations between genus two and genus one theta-functions of a similar spirit, but with different physical and mathematical details, were proved in [33][34][35]. We postpone such investigations to the future. 10 The dimension is typically linear in k (e.g. the space of modular forms on SL(2, Z) has dimension k/12 up to order one corrections), and our computations should cover these quite easily.

Summing over spin structures and the thermal entropy relation
In this section we consider the sum over all spin structures of the higher-genus result, Equation (2.16). As is well-known, there are 2 2n spin structures and one is expected to sum over all of them. There are two immediate consequences of doing so. One is that the answer manifestly satisfies Bose-Fermi equivalence, since this is also the free boson answer at R = 1. The second is that it is modular covariant, as was shown in [12]. It only remains to demonstrate that the result satisfies the thermal entropy relation, which, in its strong form, is actually a pair of relations valid respectively as z 12 → 0 and z 12 → 1. Recall that it was shown in [12] that χ t apparently cannot, in any reasonable way, be made to satisfy these relations.
For a single Majorana fermion, the sum is proportional to: where α, β range over all the 2 2n spin structures on the higher-genus surface 11 . Recall that Ω is given in terms of C k by Equation (2.9). It is convenient to parametrise C k as follows. We have already seen that C 0 = τ . As was done previously for genus 2, let us write C k = τ + 2πi C k for k = 1, 2 · · · , n − 1. Inserting this and using: As the interval size becomes small we have z 12 → 0 and C k → 0 for all k. In this limit, the second exponential tends to 1 and: which immediately implies the small-interval part of the thermal entropy relation.
On the other hand, as the interval grows large, z 12 = 1 − with → 0, we have C k → 1 π 2 sin πk n | log |. Then the second term will be exponentially damped unless the coefficient of | log | is zero. This requires: Now each β a is independently equal to 0 or 1 2 mod 1. Of the 2 n total choices, half of them have a β a = 0 (mod 1) and the other half have a β a = 1 2 (mod 1). This means that in this limit we can write: from which the z 12 → 1 limit of the thermal entropy relation follows.
Notice that in the small interval limit the replica partition function goes over to the "uncorrelated" sum over spin structures, while in the large-interval limit it goes to the "correlated" sum 12 . This was exactly the behaviour argued in [12] to satisfy the thermal entropy relation. Here the essential point is that we did not put it in by hand, rather it emerged as a property of the spin-structure-summed higher-genus Θ-function. 12 The latter statement follows by reverse applications of Equations (3.5) and (3.4) in [12] (one has to make the obvious change from the Dirac fermion case studied there to the Majorana fermion case above by taking a square root).

Conclusions and Outlook
Let us first summarise the part of our result that makes no reference to modular invariance.
Suppose one wants to calculate the n'th Rényi entropy of free fermions on a circle at finite temperature, with fixed fermion boundary conditions around the space and imaginarytime axes (the former is up to us, while the latter should be anti-periodic). Then the twist-operator method of [2] provides an answer in terms of Jacobi θ-functions, while the partition function on the genus-n replica Riemann surface provides another answer in terms of Siegel Θ-functions. We have stated a precise identity, Equation (2.24), which, if true, implies the equivalence of these two answers. We have provided some evidence for this identity for arbitrary n, and stronger evidence for n = 2. One could even turn things around and argue that free fermion theory provides the rationale, or "physics proof", of our identity 13 . Nonetheless it should be possible to work out a rigorous mathematical proof.
Now if we want to compute the nth Rényi entropy of a modular-invariant CFTfor example the Ising model-using the free fermion description, then it is clear that a sum over spin structures (fermion boundary conditions) is required [11,12]. Performing such a summation on the twist-operator computation of [2] does not provide a consistent answer compatible with physical requirements like Bose-Fermi equivalence and the thermal entropy relation. By contrast, if we sum the genus-n replica partition function over all spin structures in genus-n, we do get a consistent answer and we claim this is the correct answer for the Rényi entropy of the system. This does not exclude the possibility that some twist-operator computations, so far not performed, could generalise that of [2] to provide the complete and correct answer without recourse to the higher-genus replica partition function 14 . Of course any proposal for such a computation must agree with the highergenus replica partition function for each choice of spin structure. That would require a new set of identities generalising Equation (2.24) away from the diagonal replica spin structure.
On the way, we showed that the Rényi entropy at finite size and temperature, after removing a universal factor, transforms like a weak Jacobi form whose weight and index we obtained. In particular it is not periodic in the size of the interval or the inverse temperature. This is to be expected, since otherwise the difference between the largeinterval and small-interval entanglement would be zero, contradicting the thermal entropy relation. In fact we demonstrated a periodicity under n-fold multiples of the basic shifts.
This means that our answer contains not just the Rényi entropy but also its analytic continuation to a region where the "entangling interval" wraps the basic torus one or more 13 We thank Edward Witten for this observation. 14 We thank Matthias Gaberdiel and Shiraz Minwalla for this suggestion.
times. It would be interesting to understand the physical meaning of this more general quantity.
Finally, one may hope that our understanding of the higher-genus replica surface paves the way for the study of Rényi and entanglement entropies for other 2d conformal field theories at finite size and temperature (in this context see also [24,25] A Periodicities of χ g and χ t In this appendix we compute the periodicities of the expressions χ g (z 12 , τ ; α, β) and χ t (z 12 , τ ; α, β) defined in Eqs. (2.22) and (2.23) respectively. In fact neither side is periodic under the "simple" translations z 12 → z 12 + 1, z 12 → z 12 + τ . Rather, when n is odd we show that the two sides are perioic under an n-fold shift z 12 → z 12 + n, and quasi-periodic under the other n-fold shift z 12 → z 12 + nτ . We will find that the expressions turn out to be Jacobi forms in the variable Z = z 12 n , which is (quasi)-periodic under the standard translations Z → Z + 1, Z + τ . For even n, things are slightly different. In this case the transformations z 12 → z 12 + n, z 12 → z 12 + nτ lead to a change in spin structures for both χ g and χ t , in addition to a pre-factor in the latter case. The change in spin structures, as well as the pre-factors, are the same on both sides. For genuine (quasi)-periodicity at even n, one has to further double the shifts and consider z 12 → z 12 + 2n, z 12 + 2nτ and we find that χ g and χ t transform in the same way under these transformations.

(A.2)
For the integrals A 0k , B 0k , these shifts in z have the effect of deforming the contour of integration. We can represent the effect of this deformation as follows: The sign depends on how to take the analytic extension. However in the exponential the final result does not depend on the sign. Thus, we have: and: A 0k (z 12 + nτ ) =    A 0k (z 12 ) for k = 0 e 2πi(n−k) k n z 12 e πiτ k(n−k) (A 0k (z 12 ) − nB 0k ) for k = 0 ; B 0k (z 12 + nτ ) = e 2πi(n−k) k n z 12 e πiτ k(n−k) B 0k (z 12 ) . where B is a symmetric matrix given by: The transformation property of the Θ-function under these shifts is as follows. Consider the first case with Ω → Ω + B. Then: Using the expression for B above, one can easily show that: The second equality is obvious when α = 0, but it is easy to verify that it is also true when α = 1 2 . Now when n is odd, the last line of Equation (A.9) is even and therefore it does not modify the Θ function. Hence for odd n we have proved that: It follows that: χ g (z 12 + n, τ ; α, β) = χ g (z 12 , τ ; α, β) . (A.11) We will return to the case of even n below.