Entanglement Entropy in Generalised Quantum Lifshitz Models

We compute universal finite corrections to entanglement entropy for generalised quantum Lifshitz models in arbitrary odd spacetime dimensions. These are generalised free field theories with Lifshitz scaling symmetry, where the dynamical critical exponent $z$ equals the number of spatial dimensions $d$, and which generalise the 2+1-dimensional quantum Lifshitz model to higher dimensions. We analyse two cases: one where the spatial manifold is a $d$-dimensional sphere and the entanglement entropy is evaluated for a hemisphere, and another where a $d$-dimensional flat torus is divided into two cylinders. In both examples the finite universal terms in the entanglement entropy are scale invariant and depend on the compactification radius of the scalar field.


Introduction
Quantum entanglement refers to a correlation of a purely quantum mechanical nature between degrees of freedom in a physical system. Consider a quantum system that can be divided into two subsystems A and B, such that the Hilbert space can be written as a direct product of the Hilbert spaces of the subsystems, H = H A ⊗ H B , and take the full system to be in a state described by a density matrix ρ. The entanglement entropy of subsystem A is then defined as the von Neumann entropy of the reduced density matrix, obtained by taking a trace over the degrees of freedom in subsystem B, i.e.
with ρ A = Tr B ρ. We will take ρ to be the ground state density matrix of the full system but our results can be extended to more general states.
Entanglement entropy is a useful theoretical probe that encodes certain universal properties of field theories describing critical systems, see e.g. [1][2][3][4][5]. A well known example in this respect is the entanglement entropy of a two-dimensional conformal field theory (CFT) [6][7][8], which has a universal logarithmic term, where c is the central charge of the CFT in question, L A is the spatial size of the subsystem A, and ε is a UV-cutoff. Here we are assuming that A is connected and that its size is small compared to the full system size L A L. The logarithmic UV behaviour of the entanglement entropy tells us that the system has long-range entangled degrees of freedom (in contrast to an area-law where short-range entanglement would mainly contribute).
In recent years considerable effort has been devoted to investigating such universal terms in the entanglement entropy of CFTs in arbitrary dimensions (see [9,10] for reviews). What will be important for us is the following general UV behavior of entanglement entropy in even D-dimensional CFTs, where Σ D−2 is the area of the (D−2)-dimensional entangling surface ∂A and L A is a characteristic length associated with ∂A. Only c 0 is universal in this expression. The other coefficients depend on the regularisation scheme used. One also finds a universal sub-leading term in odd-dimensional CFTs but in this case it is finite rather than logarithmic, see e.g. [10] and references therein. The coefficient of the universal term is a function of topological and geometric invariants, such as the Euler density and Weyl invariants constructed from the entangling (hyper)-surface [11,12]. This reflects the fact that the logarithmic term in S[A] is related to the conformal anomaly of the stress-energy tensor of the corresponding CFT.
In the present work we will study entanglement entropy, including universal finite terms (i.e. of order O(1) with respect to L A ), in a family of d+1-dimensional scale invariant quantum field theories introduced in [13]. The scale symmetry is a non-relativistic Lifshitz symmetry that acts asymmetrically on the time and spatial coordinates, with a dynamical critical exponent equal to the number of spatial dimensions z = d. These theories are closely related to the well known quantum Lifshitz model (QLM), first studied in the seminal work [14]. This is a scale invariant free field theory in 2 + 1-dimensional spacetime with a z = 2 dynamical critical exponent. It is an effective field theory for certain quantum dimer models (and their universality class) in square lattices at a critical point and involves a compactified free massless scalar field [14]. Non-trivial Lifshitz scaling is achieved via a kinetic term that is asymmetric between time and space (with higher derivatives acting in the spatial directions). The higher-derivative construction can easily be extended to free scalar field theories in any number of spatial dimensions d with z = d Lifshitz scaling. In [13] such theories were dubbed generalized quantum Lifshitz models (GQLMs) and several interesting symmetry properties were revealed in correlation functions of scaling operators. The periodic identification of the scalar field did not figure in that work but turns out be important when one considers entanglement entropy in a GQLM defined on a topologically nontrivial geometry.
A key property of the QLM and GQLM theories is that the ground state wave functional is invariant under conformal transformations involving only the spatial dimensions [13,14]. The spatial conformal symmetry is a rather special feature (the corresponding critical points are called conformal quantum critical points [14]) and it manifests in the scaling properties of entanglement entropy. In essence, the symmetry allows us to map a (d+1)-dimensional Lifshitz field theory with z = d to a d-dimensional Euclidean CFT. In the d = 2 QLM the CFT is the standard free boson CFT but for d > 2 GQLM's the spatial CFT is a higherderivative generalised free field theory. Such higher-derivative CFTs have been discussed in a number of contexts, for instance in relation to higher spin theories, e.g. [15,16], as models in elastic theory e.g. [17], in a high-energy physics setting in connection with the naturalness problem e.g. [18,19], and in the context of dS/CFT duality e.g. [20,21]. These theories are not unitary but their n-point correlation functions are well defined in Euclidean spacetime and being free field theories they have no interactions that trigger instability. Entanglement entropy and its scaling properties have been extensively studied in the QLM [1,[22][23][24][25][26][27]. 1 The replica method in the QLM was first developed in [22], where it was found that for a smooth entanglement boundary, the scaling behaviour of the entanglement entropy in the QLM (and more generally for conformal quantum critical points in (2+1)-dimensions) is of the form where c 1 depends on the regulator and ∆χ is the change in the Euler characteristic (upon dividing the the system in two), which in turns depends on the topology of the system and on the entangling surface. The above behavior follows from general expectations for the free energy of a two-dimensional CFT with boundary [29]. This result is for the entanglement entropy of a non-relativistic (2+1)-dimensional theory, however its computation starts from a ground state which is a "time-independent" conformal invariant. As the time coordinate only appears as a spectator, the final result displays features of a two-dimensional CFT. This is in line with the results of [13,14], where it was shown that equal-time correlation functions of local scaling operators in the QLM and GQLM can be expressed in terms of correlation functions of a d-dimensional Euclidean CFT.
Furthermore, by choosing a smooth partition, such that we have no contribution from the logarithm (i.e. ∆χ = 0), a further sub-leading (of order one in L A ) universal term appears in the entanglement entropy for the QLM [1], that is The universal term γ QCP , where QCP stands for quantum critical point, depends both on the geometry and topology of the manifolds [23,[25][26][27]. It depends on the geometry in the sense that it includes a scaling function written in terms of aspect ratios typical of the given subsystem, and on the topology through a contribution from zero modes and non-trivial winding modes. In this sense, the entanglement entropy of the QLM is able to capture longrange (non-local) properties of the system. In particular, γ QCP was computed by various methods for a spatial manifold in the form of a cylinder in [24][25][26][27]30], for a sphere in [27], and the toroidal case was treated in [25] by means of a boundary field theory approach. The toroidal case was further investigated in [31], where analytic results were obtained for Rényi entropies. 2 Our aim is to extend the study of these finite universal terms in the entanglement entropy to generalised quantum Lifshitz models. In particular, we analyse their scaling properties in full generality, in any number of spatial dimensions d with Lifshitz exponent z = d. For technical reasons (which we explain below) we restrict attention to the case of even integer d. More concretely, we obtain universal terms in the entanglement entropy in GQLM on two classes of manifolds. On the one hand, we divide a d-dimensional sphere into two d-dimensional hemispheres, on the other hand we consider a d-dimensional flat torus, sliced into two d-dimensional cylinders. 3 Our computations are purely field theoretical. The theories we consider represent rare examples of non-relativistic critical theories, for which entanglement entropy can be obtained analytically. We view them as toy-models where we can explore quantum entanglement for different values of the dynamical critical exponent z. Further motivation comes from a puzzling aspect of Lifshitz holography, where one considers gravitational solutions that realise the Lifshitz scaling (1.4), see e.g. [34,35]. In AdS/CFT the entanglement entropy is computed by means of the Ryu-Takayanagi formula [36,37]. The usual working assumption is to apply the RT prescription also in Lifshitz holography, but in a static Lifshitz spacetime the holographic entanglement entropy does not depend on the critical exponent at all (see e.g. [38]). From a field theory point of view, however, one expects the higher-derivative terms 2 For related studies of the scaling properties of entanglement entropy for a toroidal manifold in scale invariant (2+1)-dimensional systems, see also [32,33]. 3 Here a d-dimensional cylinder refers to the product of an interval and a (d−1)-dimensional flat torus.
to dominate at short distances, and thus the UV behaviour of entanglement should reflect the value of the dynamical critical exponent. The absence of z from the holographic entanglement entropy is puzzling if Lifshitz spacetime is the gravitational dual of a strongly coupled field theory with Lifshitz symmetry. Similar considerations motivated the work in [39][40][41]. While we clearly see a dependence on z, it is difficult to compare our results directly to those of these authors as we are not working within the same class of field theories and we focus on universal sub-leading terms rather than the leading area terms.
The paper is arranged as follows. In Section 2 we briefly review the construction of generalized quantum Lifshitz models and extend their definition to two specific compact manifolds. These are higher-derivative field theories so we must ensure that the variational problem is well posed. This amounts to imposing z conditions on the variations (2.8), (2.15) and including a boundary action (2.10) for the d-torus or (2.16) for the d-sphere. In Section 3 we discuss the replica method, which we use to compute the entanglement entropy. In essence, this approach maps the problem of computing the nth-power of the reduced density matrix to a density matrix of an n times replicated field theory. The goal is to produce a result that can be analytically continued in n, in order to calculate the entanglement entropy according to (3.1). As we explain in Section 3, the replica method forces all the replicated fields to be equal at the cut, since the cut is not physical and our original field theory only has one field. There is an additional subtlety in implementing this condition due to the periodic identification of the scalar field in the GQLM and in order to ensure the correct counting of degrees of freedom we separate the replicated fields into classical and fluctuating parts. The fluctuating fields obey Dirichlet boundary conditions as well as the vanishing of the conformal Laplacian and its integer powers at the cut. Their contribution is encoded in partition functions computed via functional determinants defined on the sphere and torus respectively. The classical fields give rise to winding sectors that are encoded in the function W (n) described in Section 3. For the spherical case this contribution is simple and only amounts to a multiplicative factor √ n. For the toroidal case, the contribution from the classical fields is less trivial, and requires summing over classical vacua of the action. For higher-derivative theories some further conditions have to be implemented in the classical sector, and we argue that a compatible prescription is to use Neumann boundary conditions for these fields. At this point no freedom and/or redundancy is left, and it is straightforward to compute the sum over winding modes. We collect the contributions from the classical and fluctuating fields to the universal finite term of the entanglement entropy for a d-sphere and a d-torus in Sections 4 and 5 respectively. We conclude with some open questions in Section 6.
Most of the technical details are relegated to appendices. In Appendix A we review the computation of the functional determinant contribution for the spherical case, which was originally worked out in [42]. In Appendix B we develop an alternative expression for the formulae presented in Appendix A, which we find more transparent and better suitable for numerical evaluation. In Appendix C we compute the functional determinant contribution for the toroidal case. Finally we compute the winding sector contribution for the d-torus in Appendix D.

Generalised quantum Lifshitz models in (d+1)-dimensions
The 2+1-dimensional quantum Lifshitz model [14] can be generalised to d+1-dimensions [13]. Whenever the dynamical critical exponent z is equal to the number of spatial dimensions, the ground state wave-functional is invariant under d-dimensional conformal transformations acting in the spatial directions, extending the connection between the quantum Lifshitz model and a free conformal field theory in one less dimension to any d. We recall that in the 2+1dimensional case the scalar field is compactified [14], and below we will also compactify the scalar field in the GQLM at general d on a circle of radius R c , that is identify φ ∼ φ + 2πR c .
The ground state wave functional of the GQLM is [13] where |φ is an orthonormal basis of states in the Hilbert space made up of eigenstates of the field operator, and the partition function Z is given by We are interested in computing the sub-leading universal finite term of the entanglement entropy (1.1) in the ground state, i.e. with ρ = |ψ 0 ψ 0 |, when the manifold M is a d-sphere or a d-torus. The subsystem A will consist of field degrees of freedom on a submanifold of M.
For technical reasons we restrict attention to the case where d is an even (positive) integer.
We follow the normalisation convention of [27] and write the action as where G = det G ab (a, b = 1, . . . , d) is the determinant of the Euclidean metric on the manifold M, and P z, M is a proper conformal differential operator of degree z in d dimensions. The specific form of P z, M depends on M, as we will discuss at the end of this section. In order to have a well-defined variational problem, the action has to include a suitable boundary term S ∂M whose specific form is also given below. Note that the scalar field φ has dimension zero under Lifshitz scaling in the GQLM at general d. We find it convenient to use the shorthand g = κ 4π . 4 We note that for a flat manifold, the compactification of the field implies a global shift symmetry compatible with conformal symmetry [43,44]. This is also true for the z = d theory on the d-sphere (and more generally on any Einstein manifold) provided the action includes appropriate terms that generalise the notion of conformal coupling to a higher-derivative setting.
Let us consider how the action in (2.3) is constructed concretely for the two cases, mentioned above. To keep the discussion somewhat general, we assume that both z and d are even 4 The normalisation of the action in (2.3) and the compactification radius Rc are not independent. A rescaling of the scalar field will affect both g and Rc, while physical quantities that are independent of rescaling are expressed in terms of 2πRc √ g [43,44].
positive integers and do not insist on z = d for the time being. The case when d is an odd integer is also interesting but raises a number of technical issues that we do not address in this work. The boundary terms in the action will be important once we divide the system into subsystems and apply the replica method (cf. Section 3).
d-torus. In section 5, we consider a torus obtained as the quotient space of R d and a ddimensional lattice. The manifold is flat, and in this case the operator appearing in the action 3) is simply the z/2 power of the Laplace-Beltrami operator, For d = z = 2 this reduces to S = g (∇φ) 2 as in [14] (after integrating by parts). Varying where the partial derivatives should be understood as follows with n a an oriented unit vector normal to the boundary. We need to choose appropriate boundary conditions for the variations. One possibility is to demand that The reason behind this choice is that we will be interested in the eigenvalue problem of ∆ z/2 , for which we require the operator to be self-adjoint and to have a complete set of consistent boundary conditions. The replica method forces us to choose Dirichlet conditions on the field at the boundary (cf. Section 3) and the remaining conditions are chosen to be consistent with the self-adjointness of the operator. Equipped with (2.8), the variation of the Lagrangian reduces to (2.9) Hence, defining the following boundary action with variation given by clearly gives a well-defined variation for the total action and leads to the following equations of motion ∆ z/2 φ = 0, and ∂ 2k n δφ ∂M = 0, for k = 0, . . . , When the manifold M is a unit d-sphere, the operator in (2.3) is the so-called GJMS operator on a d-sphere [45]. In essence, GJMS operators generalise the conformal Laplacian to higher derivatives and d-dimensional curved manifolds [45] (see [46][47][48][49][50][51][52][53][54] for more references on the subject). This means that a GJMS operator of degree 2k in d-dimensions (where k is a positive integer) is constructed so that it transforms in a simple way under a Weyl transformation of the metric, G ab → e 2ω G ab , In general, the operator P 2k is well defined for k = 1, . . . d/2 for even d, in the sense that it reduces to the standard Laplacian of degree k in flat space [47]. For odd d-dimensional manifolds operators satisfying (2.13) exist for all k ≥ 1 [47,52,53].
On a unit d-sphere the GJMS operator of degree 2k explicitly reads as is the Laplace-Beltrami operator. The case of most interest to us is to consider a GJMS operator of degree 2k = d. This is known in the literature as the critical case, while k < d/2 is referred to as the subcritical case. It is straightforward to check, using (2.13), that the final action S 0 (2.3) is invariant under Weyl transformations. The factorisation in (2.14) is a general characteristic of Einstein manifolds. Since the eigenfunctions of the Laplace-Beltrami operator on a compact Riemannian manifold form an orthonormal basis, one can easily obtain the spectrum of the GJMS operator on the sphere from the factorisation above. This will play an important role later in the computation of partition functions in Section 4.
When S d is divided into hemispheres, the action S 0 [φ] has to be complemented by boundary terms. As before, in order to have a well-defined variational problem, we compute the variation of the action δS 0 , impose z boundary conditions on δφ and its derivatives, and then cancel any remaining terms against appropriate boundary terms. We impose the following boundary conditions that is Dirichlet boundary conditions on the variation δφ and vanishing of its Laplacian and its powers at the boundary. This is analogous to (2.8) for a curved manifold. The explicit where the products in the above expression are taken to be empty when the upper extreme is less than the lower extreme. Finally, the classical equation of motion is

Replica method
The entanglement entropy of subsystem A can be defined as where an analytic continuation of the index n is assumed. This definition is equivalent to the von Neumann entropy of ρ A (1.1). Following [8,55], we will use the replica approach to evaluate (3.1). At the heart of this method is to view each appearance of the density matrix ρ in Tr ρ n A as coming from an independent copy of the original theory, so that one ends up working with n replicated scalar fields. The process of taking partial traces and multiplying the replicas of ρ then induces a specific set of boundary conditions at the entanglement cuts on the replica fields.
In this section we adapt the replica trick to generalised quantum Lifshitz theories. For the QLM the replica method was reviewed in [26,27]. Our starting point is the ground state density matrix ρ = ψ 0 ψ 0 , with |ψ 0 as in (2.1). Now divide the manifold into two regions A and B and assume that the Hilbert space splits as H = H A ⊗ H B . This allows us to write the density matrix as . We then construct ρ n A by the gluing procedure represented in Figure 1. Each copy of ρ is represented by a path integral as in (3.2) with fields labelled by a replica index i = 1, . . . , n. The partial trace over field degrees of freedom with Figure 1: The gluing procedure due to the replica trick. Gluing due to the partial trace over B is represented in red, due to multiplication of the reduced density matrices ρ A in blue, and due to the final total trace in yellow. support in B gives the following reduced density matrix for the i-th replica, Multiplying together two adjacent copies of the reduced density matrix gives The δ-function coming from φ A i |φ A i+1 forces the identification φ A i = φ A i+1 , effectively gluing together the primed field from replica i and the unprimed field from replica i + 1, as indicated in Figure 1. It follows that multiplying n copies of the reduced density matrix gives rise to pairwise gluing conditions φ A i = φ A i+1 for i = 1, . . . , n − 1, and when we take the trace of the complete expression we get a gluing condition between the first and last replicas, φ A 1 = φ A n . Combining this with the gluing condition φ B i = φ B i for i = 1, . . . , n from the partial trace in (3.3), and the boundary condition where Γ denotes the entangling surface separating A and B, we see that all the replica fields are forced to agree on Γ. The final result for ρ n A is then where cut(x) is some function of the boundary coordinates. We write the denominator in (3.5) as Z n A∪B to emphasise that it contains n copies of the partition function of the original system before any subdivision into fields on A and B. In the numerator, however, the field configurations of the different replicas are integrated over independently, except that the replicated fields are subject to the boundary conditions (3.6) (up to the periodic identification φ ∼ φ + 2πR c ).
In order to take the periodic identification into account when applying boundary conditions, we separate each replicated field into a classical mode and a fluctuation, following a long tradition, see e.g. [43,44], The modes φ i,cl satisfy the following classical equations of motion and boundary conditions, are integers indicating the winding sector. The classical field determines the total field value at the entanglement cut, while the fluctuating field ϕ satisfies Dirichlet boundary conditions, In two dimensions this condition, along with the equation of motion of the classical fields, ensures that the action factorises [27], 5 The decomposition of the action is less trivial in higher dimensions but it can be achieved if the Dirichlet boundary condition on the fluctuating field at the entanglement cut is augmented by further conditions. It is straightforward to check that imposing (3.9) along with (2.8)-(2.15) on the fluctuating fields at the cut leads to a well-posed variational problem as well as self-adjointness of the operator P z,M on M. As was discussed earlier, this combination of conditions amounts to the vanishing of the Laplace operator and its integer powers acting on the fluctuating fields at the boundary. This turns out to be enough to ensure that the total action splits according to (3.10) (once again the equations of motion for the classical fields have to be used to achieve factorization). We note, that with this prescription and using the classical equations of motion, the boundary terms in the action can be written in a form that only depends on the classical part of the field, In the presence of winding modes, there remains some redundancy in the classical part of the action, as further discussed in Appendix D where we compute the contribution from the classical winding sector for the d-torus.
As a consequence of (3.10), the fluctuating modes ϕ i simply contribute as n independent fields obeying Dirichlet boundary conditions (3.9) at the entanglement cut. For the classical modes, on the other hand, we can solve for the A and B sectors simultaneously, as the boundary value problem (3.8) has a unique solution in A ∪ B Γ, up to winding numbers. At the entanglement cut only relative winding numbers matter and we can choose to write the boundary conditions for the classical fields as [26] Thus, the trace of the n-th power of the reduced density matrix reads where W (n) is the contribution coming from summing over all classical field configurations satisfying the boundary conditions (3.12). The subscript D on the integral sign is a reminder that the the fluctuating fields obey Dirichlet boundary conditions.
At this point we need to distinguish the spherical case from the toroidal one. We start by analysing the problem on the d-sphere, which turns out to be particularly simple.
d-sphere. We closely follow the treatment of the two-dimensional case in [27]. The crucial observation here is that the winding mode can be reabsorbed by the global shift symmetry of the action, i ] with constant φ 0 i , as mentioned in Section 2. Indeed, since the fields satisfying the classical equation of motion include any constant part of the total field, we can use the symmetry to rewrite their boundary conditions as We then choose φ 0 i = −2π R c ω i to cancel out all winding numbers. The boundary conditions then become φ cl i|Γ (x) = cut(x) , i = 1, . . . , n , (3.15) and the sum over classical configurations can be written as Consequently, we have for the d-sphere We can now combine the n-th fluctuating fields with support on A and B and the n-th classical field to define Notice that the effective compactification radius of Φ n is now √ nR c [26,27]. The path integral over ϕ A n and ϕ B n along with the contribution W (n) from the rescaled classical field amounts to the partition function on the whole d-sphere for the combined field Φ n , which is equal to the partition function of the original field up to a factor of √ n due to the different compactification radius, and it therefore almost exactly cancels one power of the original partition function in the denominator in (3.17), Hence, the entanglement entropy is given by 20) and the original problem has been reduced to the computation of partition functions with appropriate boundary conditions on the regions A and B and A ∪ B.
We will consider the case where the d-sphere is divided into two hemispheres. Then we only have to compute the partition function on the full sphere and a Dirichlet partition function on a hemisphere. These are in turn given by determinants of the appropriate GJMS operators. The detailed computation is described in Section 4 and Appendices A and B. d-torus. We now apply the replica method in the case of a d-torus. We cut the torus into two parts, thus introducing two boundaries which are given by two disjoint periodically identified (d − 1)-intervals (in d = 2 this is simply an S 1 ). As explained before, all the fields have to agree at the cuts Γ a (where now the index a = 1, 2 labels each cut). For the quantum fields this simply implies that they need to satisfy Dirichlet boundary conditions, that is As explained earlier, further conditions are necessary in dimensions d > 2, and we demand that the conditions (2.8) hold at the cut for the fluctuating fields.
Now consider the classical fields on the torus. We can use the global shift symmetry discussed in Section 2 to write the boundary conditions for the classical fields as with φ 0 i constant on the whole torus. As in the spherical case, we can choose the φ 0 i so that they absorb the winding numbers from one of the cuts, where ω i : = ω 2 i − ω 1 i . We are effectively left to deal with winding sectors at a single entanglement cut and since only the relative winding number between adjacent replicas matters we can eliminate one more winding number to obtain At this point, we can use the same unitary rotation U n as in [26,27], to bring the classical fields φ i,cl i = 1, . . . , n into a canonical form constructed to separate the contribution from the winding modes from the contribution from modes subject to boundary conditions given by the functions cut 1,2 (x). Concretely, we define the matrices and Hence, the sum over all the classical configurations reduces to a sum over the vector w = (ω 1 , . . . , ω n−1 ) ∈ Z n−1 and an integral over the n-th classical mode. Notice that, as for the spherical case, the n-th classical modeφ cl n has a compactification radius amplified by √ n, due to the rotation (3.27). We want to use this mode to reconstruct a full partition function on the torus, that is define where the √ n factor on the right-hand-side of (3.31) accounts for the different compactification radius. Thus, the replica method finally gives where W (n) contains the contributions from the first n − 1 classical configurations satisfying the boundary conditions in (3.29) at Γ 1(2) .
In two dimensions the classical fields are uniquely determined by the equations of motion and the boundary conditions (3.29), and thus the classical action has only one vacuum. The contribution from the winding sector is then simply given by the sum over the corresponding winding modes [26,27], (3.33) However, in higher dimensions (d > 2) the conditions (3.29) do not uniquely specify the vacua of the classical action. In other words, our construction is consistent for more than one set of boundary conditions applied on derivatives of the classical fields and the value of the boundary action depends on the boundary conditions. This is the redundancy mentioned in Section 2.
The classical field satisfies a higher-derivative equation of motion, whose general solution is parametrised by z/2 constants. The boundary condition imposed on Φ n will fix one of these constants but we need to add z/2 − 1 further boundary conditions for the classical field to fix the rest. The value of the boundary terms in the action will in general depend on the choice of boundary conditions.
In the present work we impose a generalised form of Neumann boundary conditions on derivatives of the classical fields, This prescription is compatible with the conditions imposed on the fluctuations, 35) and, at the same time, it gives a non-vanishing classical boundary action, which is important in order for the sum over winding modes to converge. The contribution from winding modes is then accounted for in any number dimentions by computing where the classical fieldsφ cl i satisfy the boundary conditions (3.29) and the Neumann conditions in (3.34). The details of the computation are reported in Appendix D.
Finally, the entanglement entropy for the d-torus is given by since the winding sector is normalised such that W (1) = 1. The computation of the partition functions for the d-torus is presented in Section 5 below.
We close this section by noting that even though winding numbers come into play across entanglement cuts in our computation, we are restricting our attention to a single topological sector of the original theory on the d-torus. Indeed, since we periodically identify the field, we could consider winding sectors on the d-torus itself, We have set m I = 0 in our calculations in the present paper but a more general study can be carried out, evaluating the contribution from winding sectors W (n, m L ) associated with an entanglement cut for each topological configuration, and then summing over the m I . The corresponding topological contributions to entanglement entropy in a scalar field theory on a two-dimensional cylinder, are obtained in [26].

Entanglement entropy on a hemisphere
In this section we calculate the universal finite terms of entanglement entropy in GQLM resulting from the division of a d-sphere into two d-hemispheres A and B by an entanglement cut at the equator as shown in Figure 2. According to the replica calculation in Section 3, we have to compute (3.20), where now A and B are the two d-hemispheres, and the bulk action contains the GJMS operator (2.14) with 2k = d. The partition function on the whole manifold (Z A∪B in (3.20)) contains a zero mode, which should be treated separately [43,44], where A d is the area of the d-sphere and the g π A d factor comes from the normalisation of the eigenfunction corresponding to the zero eigenvalue. The functional determinant det is the (regularised) product of the non-zero eigenvalues. The operator (2.14) on the unit d-hemisphere with Dirichlet boundary conditions does not have a zero eigenvalue, so the partition functions on the subsystems A and B can be directly computed via regularised functional determinants. Hence we can write (3.20) as where the D subscript on H d D indicates Dirichlet boundary conditions on the fields at the boundary of the d-hemisphere H d . At the end of the day, the entanglement entropy can only depend on the combination gA d . All factors of g/π inside functional determinants must therefore cancel out in the final result and going forward we simply leave them out of our formulas.
We now turn to the explicit computation of the functional determinants appearing in (4.2). In a series of papers [42,[56][57][58], Dowker calculates determinants of GJMS operators on spheres in any even dimension d and for any degree k ≤ d/2 via ζ-function methods. We give a selfcontained review of these calculations in Appendix A, partly to adapt them to our notation and partly to have all the results we want to use in one place. Determinants of critical GJMS operators (where the degree 2k of the operator matches d) on spheres and hemispheres are expressed in terms of multiple Γ-functions in [42]. A simplified version of these results, expressing them in terms of the more familiar Riemann ζ-function, is presented in Appendix B.
The starting point of Dowker's computation is the observation that the determinant of the GJMS operator on a d-sphere, given in terms of the spectral ζ-function, can be obtained as a sum of the corresponding determinant on a d-hemisphere with Dirichlet and Neumann boundary conditions [42,59] (again expressed in terms of spectral ζ-functions). On the hemisphere with Dirichlet boundary conditions the log-determinant of the GJMS operator is given by Here ζ is the Riemann ζ-function, and h D n , f D are given by The d k are Stirling numbers of the first kind, (z) k is a Pochhammer symbol, and M (d) is a sum of harmonic numbers and generalized Bernoulli polynomials whose explicit form is not important to us, as it cancels in the final expression for the entanglement entropy. The derivations of h D n and f D can be found in Appendix B.2, while the derivation of M (d) can be found in A.3, its explicit form is given in equation (A.58). These functions may seem quite complicated at first sight, but they all consist of well understood algebraic functions that can easily be evaluated using a computer. For the determinant of a critical GJMS operator on a hemisphere with Neumann boundary conditions we find a similar result with h N n and f N given by We note that our result in (4.6) differs from [42] by a sign in the term log(d − 1)!. This is because we treat the zero mode separately as is apparent in (4.1) and (4.2).
As mentioned above, the log-determinant on the whole sphere is the sum of the logdeterminants on the hemisphere with Dirichlet and Neumann boundary conditions [42], (4.9) With an eye towards the entropy formula (4.2), we express the ratio of determinants as where, using the properties of the binomial coefficients, one can write h n in the the following form Putting everything together, we obtain a surprisingly simple expression for the entanglement entropy of a hemisphere, with h n (d) given above in (4.11). For dimensions d = 2, 4, 6, and 8, in the critical case z = d, the entropy is given explicitly by and more values are plotted in Fig. 3. The two-dimensional case agrees with the result presented in [27]. Notice that the logarithmic term depends on the product R c √ g, which is independent of rescaling of the fields. Hence, in the case of a d-sphere cut into two dhemispheres, the finite universal terms of the entanglement entropy (4.12) are constant, they only depend on the physical compactification radius R c √ g of the target space, which appears in the above expression through the zero modes.
Explicit expressions for the subcritical case, when z < d for z and d both even integers, make apparent the relative simplicity of the results (4.12). Indeed, in sub-critical cases several additional terms involving derivatives of Riemann zeta functions appear. Tremendous simplification occurs when z = d to produce the neat expressions in (4.13). The specific expressions  Figure 3: The universal finite term (4.12) in the entanglement entropy of GQLM on a hemisphere plotted against the number of spatial dimension d (which is equal to the critical exponent z). We normalise S[H d ] with respect to the two-dimensional case, and set g = R c = 1.
for the functional determinants in subcritical cases were computed originally in [42], and are included in Appendix A.
The result in (4.12) only depends on "topological data" represented by the scale invariant compactification radius of the target space and not on other geometric features. One might object that this is because we initially set the radius of the d-sphere to one, and thus our computations are insensitive to the geometry. Indeed, as mentioned in the Introduction, for smooth entangling cuts in even-dimensional CFTs, the entanglement entropy is expected to have a universal term proportional to the logarithm of a characteristic scale of the system with a constant of proportionality which depends on the central charge and on the Euler characteristic. It can be checked that introducing a radius R of the d-sphere in our problem modifies the above results by adding a term proportional to ∆χ log R , (4.17) where ∆χ is the change in the Euler characteristic due to dividing the d-sphere along the entanglement cut. For the two-dimensional case this was understood in [1]. Just as for a two-dimensional sphere, the change in the Euler characteristic vanishes for the chosen entanglement cut (while having a non-smooth entangling surface can introduce further universal logarithmic terms). Indeed, on a non-unit sphere all eigenvalues entering our determinants are rescaled, and upon regularising this contributes, instead of equations (A.2), (A.3). Including the contribution coming from the normalisation of the zero-mode this would leave us with but it is straightforward to check, using (A.45) and (A.54), that this combination vanishes.
In fact, Dowker's construction of the determinant for the sphere as sum of determinants on hemispheres with Dirichlet and Neumann boundary conditions makes this quite transparent, since the spectral ζ-function in the Neumann case is nothing but the Dirichlet one after subtracting the zero mode. 6 Finally, we should stress that the sub-leading universal terms as (4.17) (which vanish here due to the chosen entanglement surface) are those expected in a ddimensional CFT. The quantum field theory we are considering lives on a (d+1)-dimensional manifold, and yet due to the enhanced d-dimensional symmetries in the critical d = z case, it has entanglement properties typical of d-dimensional CFTs.

Entanglement entropy on cut d-torus
We now turn our attention to the sub-leading universal terms in the entanglement entropy on a flat d-dimensional torus with circumferences L 1 , . . . , L d , The two-dimensional case is shown in Figure 4. The replica method for the entanglement entropy on the torus was discussed in Section 3, and it requires us to compute (3.37), where the winding sector contribution is given by (3.36), with the classical fields satisfying the equations of motion and boundary conditions expressed in (3.29) and (3.34). For the d-torus, the bulk and boundary terms in the action are given by (2.5) and (2.10), respectively. The operator P d,T d in (2.4) is simply an integer power of the Laplacian. We first compute the quantum contribution to the entanglement entropy arising from the partition functions in (3.37), and after that we tackle the winding sector contribution. All the detailed calculations of functional determinants are relegated to Appendix C, and those regarding the winding sector to Appendix D. In this section we collect the results and discuss some interesting limits.
The operator P d,T d has a zero mode and as result the torus partition function is given by where A d is the area of the d-torus. As was the case for the sphere, the g π factor in the determinant only amounts to a rescaling of the torus to which the entanglement entropy is not sensitive, and we can ignore it in our calculations. On the d-cylinder with Dirichlet boundary conditions, on the other hand, there is no zero mode and we can write the (sub-leading terms of) entanglement entropy as The required functional determinants are evaluated in Appendix C.
By means of equations (C.31a) and (C.6), we find that the determinant on the full torus is given by where the primed sum indicates the omission of the zero mode, K ν (z) is a modified Bessel function of the second kind and Ξ d−1 = diag 2π/L 2 2 , . . . , 2π/L d 2 is a diagonal matrix.
We have explicit expressions both for the spectral ζ-function on the torus in (C.13) and its derivative evaluated at s = 0 in (C.18), but at this stage we find it more convenient to use the above expression, and only insert explicit formulae at the end, after some cancellations.
where L = L A for Y A and L = L 1 − L A for Y B . We can rewrite the difference between the log-determinants as where the parameter u = L A /L 1 characterises the relative size of the two d-cylinders.
The explicit expression for ζ and a relabelling of the sides L i . As discussed in Appendix C, despite its appearance the above expression is rather convenient to handle, thanks to the fast convergence of the modified Bessel functions contained in the auxiliary function G. The derivative of the function G with respect to s, evaluated at s = 0, is given by where we have used the explicit expression (E.4) for the modified Bessel function K − 1 2 .
As an explicit example of the above result, the determinant ratio for z = d = 2 is explicitly given by where we used (C.20) and (C.29) and introduced the notation τ k = i L 1 L k+1 , for k = 1 , . . . , d−1, for the aspect ratios of the general d-torus.
For the winding sector, the computations are detailed in Appendix D. The end result, given in (D.13), is with Λ z given by where B z are the Bernoulli numbers. For instance, in d = 2 we have Finally, putting together the contributions from the functional determinants and the winding sector, (5.7) and (5.10) respectively, we get the following (rather long) expression for the entanglement entropy (5.3), It can be verified that the entanglement entropy is symmetric under the transformation u → 1 − u [33] as required for a pure state of the full system.
For the special case of d = z = 2 we obtain (using (5.9) and (D.16)) where The final result looks relatively simple due to some cancellations between the classical and quantum contributions. To our knowledge, this is the first time that universal finite terms in the entanglement entropy on a torus have been obtained in closed form using path integral methods, even for the two-dimensional case. They have been computed numerically in [30] and by means of a boundary field method in [25].
We will now check some interesting limits of our general expressions.
Halved d-torus. The first simplifying special case that that we consider is when the torus is divided into two equal parts: The contribution from the functional determinants (5.7) simplifies tremendously, leaving only a single term, 16) and we obtain for the universal terms in the entanglement entropy where now In d = 2 this reduces to since the contributions from the Dedekind-eta functions in (5.14) cancel against each other when u = 1/2. Moreover, Λ 2 is given here by Thin d-torus. In d = 2, the infinitely thin torus limit (sometimes called also the long torus limit) amounts to |τ 1 | 1 and u fixed. It can be helpful to think of this limit as L 2 → 0 while all the other lengths (L 1 , L A ) are kept fixed. In this case, the contribution from the integral in the expression (5.14) is exponentially suppressed, and moreover, we have the asymptotic behaviour (E.7) for the Dedekind η function. It is then clear that all the contributions from the Dedekind η-functions vanish in this limit, and we are left with the simple result which agrees with [25]. In this limit the entanglement entropy for the thin torus is twice the entanglement entropy for the thin cylinder [25], since the entropy still carries information about the two boundaries of the torus.
We can take a look at the same limit for the d-torus. In the d-dimensional case we assume that L 1 , L A are fixed and of order one, while all the other sides are approaching zero, that is L 2 , . . . , L d → 0. There is an ambiguity in how to take this limit, so as a first step we consider the case which can also be written as for even d. We can rewrite the term more elegantly as a function of the aspect ratios, 7 Finally, the integral over k in (5.13) is also exponentially suppressed and will not contribute to the final expression. Then, keeping only the most divergent term according to (5.22), we obtain Similar limits were discussed in [33] for the Renyi entropies of 3+1-dimensional relativistic fields theories with various twisted boundary conditions. Except for having the same powerlaw divergence, our results appear not to agree with their findings. The comparison is tricky though, as there are effectively three length scales in the d = 4, and since we are looking at the regularised entanglement entropy we do not have an explicit cut-off as in [33].
We should stress that when d = 2 all the sums in ζ The thin sliced d-torus. In d = 2 this limit corresponds to L A → 0 while all the other length scales involved remain fixed, that is u → 0 while |τ 1 | is kept fixed. The integral in (5.14) can then be evaluated, for instance by means of the Poisson summation formula (E.2), and at leading order it gives 1 2 log u. Considering only the leading term in the expansion of the Dedekind function (E.6), we obtain, for u → 0, which agrees with the entanglement entropy for the infinite long and thin sliced cylinder computed in [27]. Indeed, in this limit the torus and the cylinder are indistinguishable at leading order.
We can proceed with similar arguments in higher dimensions, assuming u → 0 while all the aspect ratios |τ i |, with i = 1, . . . , d − 1 are kept fixed. In order to simplify the computation we assume all the aspect ratios to be equal, |τ i | = σ for all i = 1, . . . , d − 1. Then, the leading divergent terms are contained in G (0, 2L A , L 2 , . . . , L d ) in (5.13), and by estimating the d−1-dimensional sum in G (0, 2L A , L 2 , . . . , L d ) (cf. (5.8)) with an integral we obtain the following leading behaviour for the entanglement entropy, where κ d is a numerical coefficient that depends on the number of dimensions d. Similar behaviour was obtained for the three-dimensional torus in conformal field theories in [62] (see also [63]), and also in [64] from a holographic approach.
The wide d-torus. As our final example, we consider the so-called wide torus limit, that is when the directions transverse to the cut are very large while L A , L 1 are kept fixed. Let us start by considering this limit for d = 2. This means that |τ 1 | → 0 while u is kept fixed.
Using the expansion of the Dedekind-eta function (E.6), we see that the term containing the logarithm of the ratio of Dedekind-eta functions in (5.14) produces the leading divergence. Hence, from the general expression for the entanglement entropy in d = 2 (5.14), we obtain This asymptotic behaviour is also expected for the universal function of the Renyi entropies of the two-dimensional torus, cf. [33] and references therein, and was found in holographic CFTs in [32].
In higher dimensions we can consider the limit when u is kept fixed, and all the transverse directions are very large compared to L 1 , L A , but all the aspect ratios approach zero at the same rate, that is |τ i | = ε, with i = 1, . . . , d−1 and ε → 0. In this case, the expressions in (5.13) simplify, and, as in the two-dimensional case, the leading divergent contribution is contained in the functions G (0, L, L 2 , . . . , L d ) (cf. (5.8)), where L can be L 1 , 2L A or 2(L 1 − L A ). Using a similar expansion as performed in the thin sliced torus limit, we obtain where f d (u) is a function symmetric under the exchange u → 1 − u.
In the above discussion, the case u = 1 2 is special for any dimension d, since the function f d (u = 1/2) = 0, so that the sub-leading but still diverging terms become important. Looking directly at (5.17) and (5.19), there is no contribution now coming from G (0, L, L 2 , . . . , L d ), and the next divergent term is logarithmic in the aspect ratios τ i , which in the two-dimensional torus is entirely coming from the integral in (5.19), while in higher dimensions it receives contributions also from the area term log √ A d and ζ T d−1 (0). In our simplified limit where all the ratios |τ i | approach zero at the same rate, we see that This is consistent with the findings of [33], where for the z = 2 free boson field theory in 3+1 dimensions the universal function of Renyi entropies J n satisfies the relation lim |τ |→0 J n (u = 1/2, |τ |)|τ | = 0. 8 This clearly holds in our case since the universal term of the entanglement entropy has a logarithmic divergence. Rather different behaviour was observed in free twodimensional CFTs for Renyi entropies [33] and also in holographic CFTs in two and three space dimensions for entanglement entropy [64], where also for u = 1/2 the universal part continues to have a power-law divergence similar to (5.28) and (5.29). The disagreement was already observed in [33].

Discussion
In this work we have analytically computed the universal finite corrections to the entanglement entropy for GQLMs in arbitrary d+1 dimensions, for even integer d, and on either a d-sphere cut into two d-dimensional hemispheres or a d-torus cut into two d-dimensional cylinders. GQLMs are free field theories where the Lifshitz exponent is equal to the number of spatial dimensions, and they are described in terms of compactified massless scalars. When d = z = 2 the GQLM reduces to the quantum Lifshitz model [14], and our findings confirm the known results of [25][26][27]. The calculations are performed by means of the replica method. Caution is required when performing the cut as the massless scalar field in the GQLM is compactified. It is useful to discern between the role of the fluctuating fields and the classical modes. In essence, the periodical identification mixes with the boundary conditions imposed at the cut on the replicated fields [23,[25][26][27]30], and disentangling the winding modes from the rest leads to an additional universal sub-leading contribution to the entanglement entropy. The fluctuating fields satisfy Dirichlet conditions at the cut (as well as further conditions imposed on their even-power derivatives), while the classical fields take care of the periodic identification. The contribution from the fluctuating fields comes from the ratio of functional determinants of the relevant operators, which we compute via spectral ζ-function methods. The classical fields contribute via the zero-mode and via the winding sector summarised in the function W (n).
For the spherical case, the full analytic expression of the universal finite terms turned out to be a constant, depending only on the scale invariant compactification radius, cf. (4.12). This is a consequence of the presence of zero modes, and their normalisation.
For the toroidal case, the story is rather rich. The general expression is (5.13), while (5.17) is valid when we cut the torus by half. In both cases, the universal term is comprised of a scaling function, which depends on the relevant aspect ratios of the subsystems, and a constant term, which contains the "physical" compactification radius. The last one comes from the zero mode of the partition function of the d-torus as well as from the winding sector. We considered various limits, such as the thin torus limit, which results in the simple expressions (5.21) and (5.25) in two and d dimensions respectively, the thin sliced d-torus, cf. (5.26) and (5.27), and finally we examined the wide torus limit, cf. (5.28), (5.29), and (5.30) where the last expression is valid for u = 1/2. Notice that in the toroidal case, where the winding sector is non trivial, it also contributes to the scaling function. For example, in the thin torus limit its contribution is decisive in order to cancel divergences and leave a finite result, cf. (5.21) and (5.25). Our findings confirm expectations from the study of the (2+1)dimensional QLM [24][25][26][27], that also for critical non-relativistic theories entanglement entropy can encode both local and non-local information of the whole system. Moreover, our results give substance to the field theoretic intuition that entanglement entropy should depend on the dynamical critical exponent, see also [39][40][41] for analogous results in this direction. The specific dependence is rather non-trivial already in the simplest spherical case. The next step would be to extend our analysis to even-dimensional spacetime, that is when d is an odd integer. Progress has been recently made in e.g. [65] (and references therein), in computing determinants of GJMS operators on odd-dimensional spheres, see also [66] for a different approach based on heat kernel techniques. It would also be interesting to broaden our study of GQLMs by examining entanglement entropy for non-smooth entangling surfaces. In the QLM for non-smooth boundaries, sharp corners source a universal logarithmic contribution (i.e. log L A ), with a coefficient that depends on the central charge and the geometry of the surgery [22,29]. For d-dimensional CFTs and in presence of non-smooth entangling surfaces further UV divergences also appear whose coefficients are controlled by the opening angle. In particular, for a conical singularity, the nearly smooth expansion of the universal corner term (seen as a function of the opening angle θ) is simply proportional to the central charge of the given CFT [67][68][69][70][71][72][73][74]. 9 Another relevant direction to pursue is the study of post-quench time evolution of entanglement in these systems, see e.g. [77] for recent results in QLM and [78] for Lifshitz-type scalar theories. 10 In this respect, GQLMs can provide an interesting and rich playground where one can answer many questions in a full analytic manner.

A The determinant of GJMS operators on spheres and hemispheres
In this appendix we review the calculation of the determinants of GJMS operators on spherical domains, originally performed by Dowker in [42], and rewrite his results in a way we find more transparent. Building on his previous work, in particular [56] and [57], Dowker writes the following expression for the spectral ζ-function of a GJMS operator of degree 2k, which we denote by P 2k , on the hemisphere, Our task therefore is to evaluate Z d (0, a D/N , k). In order to do so, we first review the definition and key properties of the Barnes ζ-function. We then turn to the evaluation of the spectral ζ-function of one of the factors in (A.1) and finally put everything together to get the desired log determinants.

A.1 Definition and properties of the Barnes ζ-function
The Barnes ζ function is defined for s > d as For the special case of d = 1 = (1, . . . , 1) we use the simplified notation, l (a) are generalized Bernoulli polynomials. For future convenience we define such that for k = 1, . . . , d we can simply write Res s=k ζ d (s, a) = R d (k, a). In particular this means that we have the following expansion around s = 0 for k = 1, . . . , d, Calculating ∂ s sζ d (k + s, a)| s=0 using the analytic continuation (A.6) one finds that where the derivative of R should be read as R d (k, a) = lim s→0 ∂ s R d (k + s, a), and ψ is the digamma function. For non-negative k, ζ d (−k, a) can be evaluated using its analytic continuation (A.6). On page 151 of [80] the following expression is given where again B Below, we derive the following two identities

A.2.1 Dirichlet boundary conditions
We now proceed with the derivation of (A.15). Since we will only be dealing with Dirichlet boundary conditions throughout this section, we will omit the D-index in a D and only reintroduce it when we reach the final result. The first step in the evaluation is to perform a binomial expansion of ζ D (s). The expansion converges because (a + m · d) 2 ≥ α 2 for d = 1 in both the Dirichlet and Neumann cases, with equality only arising for the zero mode of the Neumann case, which is omitted anyway. We thus get 2s + 2r, a), where we used the definition of the Barnes ζ-function for d = 1 given in (A.5) and rewrote the Pochhammer symbol for r > 0 as (s) r = s(s + 1) · · · (s + r − 1) = : sf (s) to make the behaviour around s = 0 more transparent. We also note that f (0) = (r − 1)!. Considering that ζ d (2r) has simple poles at r = 1, . . . , d/2 and converges for higher r, we can immediately write down the following expression at s = 0 The next step is to evaluate ζ D (0). In order to do this, we first take a partial derivative with respect to s at general values of s and then let s → 0, We note that f (0) = (r − 1)!H r−1 , where H r−1 is a harmonic number, and remind the reader that for r = 1, . . . , d/2 and small s and similarly Thus at s = 0 we can write Putting these expressions together we find There is no contribution to ζ d after r = d/2, since all such terms in the sum get set to 0 by the s factor.
The next task is to compute the remaining infinite series in (A.23). In order to do this we first note that ζ d (s, a) admits the following integral representation for s > d Using this we can rewrite where we introduced a convergence parameter σ in order for the integral representation to make sense. For the σ → 0 limit we will need to use the analytic continuation of ζ d . We now take a closer look at the different parts of the expression above around σ = 0. The ζ d functions on the first line of (A.25) are convergent at σ = 0 while the Gamma function has a simple pole. Expanding order by order in σ gives Carrying out the corresponding expansion on the second line of (A.25) gives To get a finite end result we need the pole to vanish when we add (A.26) and (A.27). 11 This is equivalent to the condition 2r, a).
We can then take the limit and write 2r, a), where H 2r−1 ≡ ψ(2r) + γ is a harmonic number. For clarity and later convenience we write down the following identity resulting from the above discussion in the case of Dirichlet boundary conditions explicitly Putting everything back together into (A.23) we can write for the complete ζ-function in the Dirichlet case where H 1 (r) is a special case of defined for any integer n ≥ 1. Higher values of n will be relevant when carrying out the corresponding analysis for a zeta function of the form ζ d (2ns + 2r, a) instead of ζ d (2s + 2r, a).

A.2.2 Neumann boundary conditions
Calculating the Neumann case is equivalent to setting a N − α = ε and letting ε go to 0. In this limit Barnes [81] calculated that is called a Γ-modular form and obeys the following identity involving the multiple Gamma function .
The multiple Gamma function Γ d (a|d) is defined, and we write Γ d (a) ≡ Γ d (a|1). Equipped with these tools we can compute the derivative of the ζ-function in the Neumann case (A.14), where ζ N is just ζ D but with a D replaced by a N . We can also write the following Neumann version of the identity (A.30) 2r, a N ). (A.37)

A.3 Determinant of GJMS operators
We are now prepared to tackle the calculation of the log-determinants. In section A.3.2, we discuss the Neumann case, which is a little more involved, as the zero mode makes the log-determinant ill-defined in the critical case k = d/2. For the subcritical case k = 1, . . . , d/2 − 1 we get while for the critical case we get the slightly more complicated result where M (d) = M (d, d/2) and ρ d+1 is a Γ-modular form as defined in (A.34).
In section A.3.3 we assemble these results to construct the determinant of GJMS operators on a sphere, summarised in (A.60).

A.3.1 Determinant on the hemisphere with Dirichlet boundary conditions
We will now generalise the results of the last section. Our starting point is the spectral ζfunction defined in (A.1). As before, we start by doing a binomial expansion for every term in the product, r=(0,...,0,r i ,0,...,0) We remind the reader that α j = j + 1 2 , d = 1, and a D = (d + 1)/2 for Dirichlet boundary conditions. For simplicity we write a instead of a D throughout this section.
The remaining sums, denoted by R(s) above, will vanish for Z d and Z d around s = 0. This is easy to see, as (s) r = O(s) while the ζ d functions will contribute with simple poles and ζ with double poles. When more than two Pochhammer symbols are present, they overcome the poles and give zero in the limit. We can make further simplifications in some of the sums above by noting that the r i are dummy indices and rewrite (A.41) as  2r, a) .
Using equation (A.28) we can rewrite this as By virtue of (A.12) evaluated at k = 0, this can in turn be written in terms of generalised Bernoulli polynomials [42], Evaluating the derivative at s = 0 is also relatively painless now. For the first two terms the result follows directly from the discussion in the previous sections, in particular from equation (A.30) in the Dirichlet and (A.37) in the Neumann case. In the Dirichlet case we thus have d, a, k), where we defined M 1 in the last line to simplify the notation. To calculate the contribution from the double sum in (A.42) we remind the reader that ζ d (s + n, a) = 1 s R d (n, a) + O(s 0 ) and ζ d (s + n, a) = −1 s 2 R d (n, a) + O(s) for n = 1, . . . , d. Since (s) r = : sf (s) = O(s) with f (0) = (r − 1)!, this means that the only terms in the sum that will survive are those for which 2r + 2r ≤ d. Also noting that ∂ s (s) r s=0 = (r − 1)! we get the following result  n (α−x) [80] of generalized Bernoulli polynomials implies, according to (A.9), that R d (2r, a N ) = R d (2r, a D ) for even dimension, which in turn implies that M (d, k) does not depend on whether we chose a D or a N . It is possible to further simplify the result by using (A.34) and noting that α j +n = α j+n . This induces a telescope-like cancellation in the product, Putting everything together we get the following expression for the log-determinant which is valid for k = 1, . . . , d/2, i.e. in the critical as well as subcritical case.

A.3.2 Determinant on the hemisphere with Neumann boundary conditions
In the subcritical case there is no zero mode and we can continue from (A.50), inserting the appropriate a for Neumann boundary conditions, i.e. a N = (d − 1)/2. Then the expression for the spectral ζ-function is as in (A.45), thanks to the above mentioned property of the generalised Bernoulli polynomials [42], and the subcritical expression for the determinant is In the critical case k = d/2 the calculation is more subtle. In order to get a well-defined expression, we must subtract the zero mode contribution from the sum, We can now use (A.48) as well as (A.33), to write down the following expression, where we denote ε = a N − α k−1) . The logarithmic divergence cancels and we obtain where we used the fact that a N + α k−1 = d − 1 and k−2 j=0 log a 2 N − α 2 j = log (d − 2)! for Neumann boundary conditions. As before, we can use identities involving the multiple Gamma function to induce cancellations. This leads us to the final result in the critical case, where we used ρ d = Γ d+1 (1), and the abbreviated notation M (d, d/2) = : M (d). Since the critical case is the most relevant to this paper, we give the explicit form of M (d), (A.58) As mentioned in the main body, our result differs by Dowker's result due to the different sign in front of the log(d − 1)! term because we are considering the functional determinant of the GJMS operators on the sphere with the zero mode removed.

A.3.3 Determinant on the Sphere
The log-determinant on the sphere is obtained by adding the log-determinants on the hemisphere for Dirichlet and Neumann boundary conditions. In the critical case, the spectral ζ-function for the sphere at s = 0 is given by (A.45) and (A.54), [42], so that 59) and the log-determinant on the sphere is thus given by (A. 60) In the subcritical case we instead have and the functional determinant reads (A.62)

B An alternative form for determinants of GJMS operators on spheres and hemispheres
In this appendix we will show how to rewrite Dowker's expressions (A.51) and (A.57) for the log-determinants of GJMS operators on hemispheres as where ζ is the Riemann ζ-function, and h n and f are functions that we derive and that depend on the boundary conditions, the dimension, and the degree of the GJMS-operator. The spherical case then follows directly as the sum of Dirichlet and Neumann hemispherical results,

B.1 Rewriting ζ d in terms of the Riemann ζ-function
In [82], Adamchik gives a closed form of the Barnes ζ-function in terms of a series of Riemann ζ-functions. His calculation is summarized by equations (14), (17), and (23) in the reference, which in our slightly different notation and after a bit of rearranging read where we note that the G d (z) are multiple Gamma functions with a different normalization compared to the one used by Dowker. We have the following explicit closed forms for all the parts of the above expression, For our purposes, it is not necessary to know how the polynomials P k,d (z) are defined, it suffices to know that they satisfy [82] d−1 We are only interested in the special case of the above formula, for which z is a positive integer. For z = 1 it is clear that log We can now define the following quantity It is further clear from ζ(s) = ζ(s, 1) that A(d, 1) = 0. Another very useful simplification comes from rewriting the R d−k sum in (B.3) as where we define D k (d, z) in the last line as With these definitions it possible to write down the ζ d (0, z) in the following very simple way

B.2 Determinants in terms of Riemann ζ-functions
Dirichlet boundary conditions on the hemisphere We can now use the above technology to rewrite the log-determinant (A.51) on the hemisphere for Dirichlet boundary conditions, (B.12) Here we have defined two functions (B.14) Since the critical case k = d/2 is of particular interest to us, we give the explicit form of h D n and f D in that case, with the functions h N n (d, k) and f N (d, k) defined in the last line as The critical case (A.57) is slightly harder, as we also need to rewrite log ρ d+1 . In order to do this, we use the following formula from [82], Reminding the reader that d d+1 = 0 we can thus write

C Functional determinants on the flat d-torus
In this appendix we calculate functional determinants of the Laplacian on the flat torus (Appendix C.1) and on a cylinder obtained by cutting the torus along a cycle and imposing Dirichlet boundary conditions along the cut (Appendix C.2). After that, in Appendix C.3, we derive the functional determinant of k-powers of the Laplacian both on the flat torus as well as on the cylinder.

C.1 Determinant of the Laplacian on the flat torus
The d-dimensional flat torus can be defined as the hyperinterval with side lengths given by L i and opposing sides identified, Since the manifold is flat, the Laplace-Beltrami operator is just the standard Laplacian on Euclidean space.
We compute the functional determinant for the Laplacian on the d-torus via a spectral ζ-function method. The eigenvalues of the Laplace operator −∂ a ∂ a , are given by As usual, the zero mode, with n 1 = n 2 = · · · = n d = 0, needs to be removed in the computation of the spectral ζ-function. In order to calculate the determinant of the Laplacian, we need to evaluate where n d : = (n 1 , . . . , n d ) is a d-vector, Ξ : = diag 2π/L 1 2 , . . . , 2π/L d 2 a d × d matrix, and the prime on the sum denotes the omission of the zero mode.
We can directly evaluate the log-determinant on the flat torus, as the analytic continuation (C.4) has only a pole at s = d 2 in the whole complex plane, This expression is recursive, and thus not in a closed form yet, however, it turns out to be useful in the evaluation of the functional determinant contribution to the entanglement entropy in Section 5.
We now solve the recursion directly and provide a closed form for the ζ-function. Our calculation of the closed form is slightly more direct than the one carried out in section 4.2.3 of [84] but the end result is the same. In order to make the calculation more transparent, let us for a while rewrite (C.4) as where the functions above summarize the information in the recursion: (C.8d) In this notation it is not difficult to see, that after k steps we get The anchor for the recursion is We note that the products of f k functions in (C.9) simplify to hence, performing the recursion for k = d − 1 steps and reinserting the definitions then gives us the following expression with the understanding that for j = 1 the product L 1 . . . L j−1 is simply 1. If we in addition insert the definition of G we can rewrite the complete expression as so the only term that survives is the one where the derivative hits the Γ-function, and we effectively set s = 0 everywhere and the Γ-function to 1. For the first sum the situation is a bit trickier, as there are two Γ-functions involved, whose poles cancel only for even j, meaning that we have to separate even and odd j for the calculation. Let us first take a look at the even j part of the first sum, that is j = 2 , we have d ds When j is odd, that is j = 2 − 1, we get d ds Before putting everything together, we note that the sums over the modified Bessel functions converge exponentially at s = 0 [83]. We can thus introduce the following notation for their limits The expression above gives us the functional determinant for a d-dimensional torus, and despite its intimidating appearance, it is rather straightforward to handle, since the modified Bessel functions hidden in S converge rapidly to zero as the integers n i increase.
It is instructive to evaluate the d = 2 case. First of all, we note that the sum on the first line of (C.18) is empty for d = 2. The remaining two sums in (C.18) consist only of the j = 1 term, giving us Using more standard conventions, see e.g. [43], introducing the modular parameter τ = iL 1 /L 2 , as well as defining q := e 2πiτ , and then performing the sum over n 1 , we obtain where η is the Dedekind η-function, defined as in (E.5). If we denote the area of the 2-torus by A, where A = L 1 L 2 in our convention, 12 then taking into account the contribution from the zero-mode, the log-determinant and the partition function on the torus can be written as which is a well known result, see e.g. [43].

C.2 Determinant of the Laplacian on the cut d-torus
We now consider cutting the torus T d L 1 ,...,L d at x 1 = 0 and x 1 = L < L 1 , as shown in figure 4. As discussed in the main body, this gives rise to two subsystems, each of which is a cylinder represented by an interval times a d−1-dimensional torus. Imposing Dirichlet boundary conditions (3.21) the eigenvalues of the Laplacian −∂ a ∂ a are µ m,n 2 ,...,n d = mπ L 2 + λ n 2 ,...,n d , m ∈ N + , n 2 , . . . , n d ∈ Z , (C. 21) with λ n 2 ,...,n d the eigenvalue on T d−1 L 2 ,...,L d . Notice that there is no zero mode now. The spectral ζ-function on this geometry thus takes the form where we schematically write λ for the (d−1)-dimensional toroidal part of the eigenvalue, and we explicitly separate the (d−1)-dimensional toroidal zero mode from the rest in the last passage. We can now evaluate the primed sum by means of the identities (E.1), (E.2), and (E.3) collected in appendix E, and we obtain Notice that λ is nothing but n d−1 Ξ d−1 n d−1 as defined in Appendix C.1. Hence, adopting the same notation here, the above term containing the modified Bessel function can be written in terms of the function G as defined in (C.5), and we can finally write the expression for the ζ-function on the cut torus as follows, (C.25) Our expression (C.25) agrees with the results of [85] obtained by contour integration.
Let us check that for d = 2 we obtain the well-known result for the log determinant of the Laplacian on a cylinder. In this case T 1 L 2 is nothing but a circle of length L 2 , and according to (C.10), we have Hence, we obtain This is a well-known result in literature, see e.g. [43,44]. It is convenient to leave α general, so that we can easily use the above results for the functional determinants for arbitrary cuts.

C.3 Determinant of powers of the Laplacian on the torus
When calculating the entanglement entropy of the GQLM on a d-torus with d even, the determinants that arise are those of even powers of the Laplacian. On the flat d-torus geometry the higher-derivative conformal operator P z is indeed just the z/2-th power of the standard Laplacian, cf. equation (2.4) in Section 2. In order to generalise our previous result, we first make the observation that, since the flat torus as well as the cut torus are compact manifolds, the spectrum of ∆ k is just given by the set of λ k , where λ is an eigenvalue of ∆ as in (C.2) or (C.21) in the case of the d-torus or the cut d-torus respectively. In particular, the spectral ζ-function corresponding to ∆ k is given by where we write schematically T for either T d L 1 ,...,L d or [0, L] × T d−1 L 2 ,...,L d , and λ for the corresponding eigenvalues (C.2) or (C.21) respectively. This leads to the simple result for the determinants [86] log det ∆ k where the critical case is found by setting k = d/2.

D The winding sector for the d-torus
The goal of this appendix is to compute the winding sector contribution (3.36) that originates from cutting the d-torus, as discussed in Section 3. In order to do so, we first need to solve the classical equations of motions for the n − 1 classical fields, (2.12), obeying the boundary conditions (3.29) as well as (3.34). As discussed in Section 3, the n-th classical field is reabsorbed into the constrained partition functions to create a free one, cf (3.31), thus here we are only concerned with n − 1 classical fields.
To facilitate reading, we list here again equations of motion and explicit conditions which the n − 1 classical fields have to fulfill. The equations of motion are given by The standard way of solving such a partial differential equation is by separation of variables, and here it suffices to separate only the first variable x 1 from the remaining orthogonal d − 1 directions y, as e.g.φ Notice that the last condition above has to be imposed either at x 1 = L A or at x 1 = −L B , or, in other words, we solve for classical fields in the region A and in the region B separately and then we glue the solutions at the boundary. As we discussed in Section 3, these conditions (D.5b)-(D.5c) are not sufficient to specify a solution of the equation of motion (D.5a), which is in general a polynomial of degree z − 1, expressed in terms of z coefficients. A choice of supplementary boundary conditions are given by (3.34), which now imply ∂ n ∆ kφcl i Γ 1 = 0 k = 0, . . . , (D.6b) We then solve the equations respectively in A and B, and glue the solutions at the two cuts, that is the whole solution is given by While the explicit value of the coefficients is quite complicated, their i dependence is simple. For instance, for the first few values of z, the functions f i,A read Notice that the functions f i are continuous everywhere on A ∪ B, but they are not differentiable at the cuts Γ 1 , Γ 2 , even though the left and right derivatives exist and are finite. This is not unusual, and it is true already in d = 2, see for example the cylindric case in [26,27]. Now we are ready to evaluate the winding sector contribution (3.36), that is the function The action is given by the bulk term S 0 (2.5) and the boundary part S ∂ (2.10), and it is clear that only the latter will contribute to W (n). As it turns out, evaluating S ∂ [φ cl i ] results in the following rather simple closed expression (notice that only the highest derivative term contributes to the action (2.10)), 13 where the B z are Bernoulli numbers. Notice that only the case z = d gives a scale invariant winding sector contribution as expected from a critical theory. With this result, we can use the analytic continuation found in appendix F of [27] to write W (n) = ω∈Z n−1 exp −π Λ z (L A , L B , L 2 , . . . , L d ) ω T · T n−1 ω (D.10) The derivative with respect to n at 1 is finally given by (D. 13) In order to make comparison in a more transparent way it is useful to rewrite Λ z in terms of the aspect ratios of the d-torus. Introducing the dimensionless parameters τ , as well as u the parameter controlling the cut, as follows 14) and noticing that L B = L 1 − L A , we can write In particular, for d = 2 and z = 2, we obtain (D.16) 13 As we discussed the classical fields are not differentiable at the cuts, however left and right derivatives exist and they are finite, so here the integrals are evaluated using left and right limits, exactly as in d = 2 dimensions.

E Useful formulae
Here we collect some useful formulae and definitions of special functions used in the paper.
The integral representation of the Gamma function immediately leads to (E.1) The Poisson summation formula is given by The integral representation of a modified Bessel function is and the explicit expression for the special case ν = − 1 2 is The Dedekind function is defined as follows η(τ ) := q 1/24 ∞